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Scope  of  Project 

To  conduct  a  basic  theoretical  research  program  in  stress  waves 
and  penetration  mechanics,  with  particular  emphasis  on  armor  plate.  To 
assist  the  experimental  program  in  this  field  being  conducted  at  the 
Arsenals  by  paraLlel  theoretical  investigations. 

Project  duration  -  from  1  January  19 63  to  15  May  1966. 

The  technical  work  of  this  project  has  been  reported  by  presen¬ 
tation  at  professional  meetings,  through  publication  in  scientific 
journals,  supported  by  reprints  submitted  to  AROD,  as  far  as  is  avail¬ 
able,  and  by  Interim  Technical  Reports.  While  recognizing  the  greater 
usefulness  of  journal  articles,  the  purpose  of  the  Interim  Reports  has 
been  to  make  the  results  available  much  sooner  than  reprints  of  publica¬ 
tions.  At  times  the  lag  has  been  reduced  by  as  much  as  a  year  thereby. 

A  second  purpose  of  the  technical  reports  has  been  to  provide  a  much 
fuller  presentation  than  is  possible  in  a  paper  and  as  a  repository  for 
data  and  details . 

The  contents  of  this  final  report  include: 

1)  list  of  output  documents, 

2)  general  summary  and  review  of  project  activities, 

5)  progress  of  last  6  month  period,  and 

4)  technical  work  of  significance  not  previously  reported. 
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Brief^Revl^j^^^Pnjject 
1.  Oblique  Impact 

One  of  the  phases  of  the  research  on  this  project  has  been 
the  analytical  study  of  tran  lent  stresses  in  a  plate  due  to  step  im¬ 
pacts  at  a  point  on  the  surface.  The  point  impact  is  a  fair  idealiza¬ 
tion  representing  a  pointed  projectile,  and  is  applicable  to  other  types 
of  disturbances  as  well.  Considerable  work  on  this  problem  was  carried 
out  on  an  earlier  AROD  project  for  normal  and  shear  impacts,  using  a 
mathematical  procedure  which  has  come  to  be  known  as  "Caignard’s  Method", 
used  extensively  in  geophysical  analyses.  The  present  effort  continued 
the  study  and  brought  it  to  conclusion  by  working  out  the  principal 
stresses  induced  for  impacts  *»t  j0°  and  60°  angles  of  incidence.  The 
results  are  reported  in  paper  no.  1  md  report  no.  1,  and  a  computer 
program  is  available  for  any  additional  numerical  data  required. 

This  method  of  analysis  is  valuable  in  that  it  follows  both  the 
dilatation  and  distortion  wavefronts  or  pulses  across  the  plates  as  well 
as  their  interaction  due  to  reflections  across  the  back  face.  This  pro¬ 
vides  a  direct  understanding,  unlike  other  methods  based  on  "normal 
modes"  or  vibrational  considerations.  Further,  it  is  able  to  assess 
realistically  magnitudes  of  the  inuuced  tensile  stresses,  which  are 
very  sensitive  to  the  reflected  waves.  For  this  reason,  analyser,  based 
on  seal-infinite  media  or  on  normal  modes  lead  to  consi  ierable  errors 
on  this  estimation.  On  the  other  hand,  Caignard’s  approach  is  not 
realistic  for  studying  actual  penetration  dynamics  (for  which  it  was 
not  intended)  and  is  limited  to  brittle- type  materials  which  are 
elastic  until  failure. 
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Since  penetration  dynamics  was  considered  the  main  theme  of 
this  project,  the  following  farther  phases  were  studied. 

2.  Plug  Formation 

This  phase  was  concerned  with  an  area  of  penetration  dynamics 
referred  to  as  "plug  formation".  This  is  part  of  the  more  general  area 
of  plastic  impact,  where  the  impacting  materials  undergo  a  type  of 
viscous  flow.  It  may  be  considered  as  intermediate  in  impact  velocity 
range  between  the  elastic  and  hydrodynamic  deformation  theories,  and 
is  characterized  by  a  paucity  of  results  on  account  of  the  uncertainties 
in  the  constants  of  the  material  and  complexity  of  the  analyses  hereto¬ 
fore  presented  in  the  literature.  A  reasonable  assumption  made  by 
Pytel  on  a  previously  sponsored  AROD  project  (also  by  some  other  litera¬ 
ture)  has  been  that  of  linear  viscosity,  so  that  the  methods  of  the 
equation  of  diffusion  or  heat  conduction  were  applicable.  Although  it 
becomes  possible  to  solve  this  equation  explicitly  with  the  given 
boundary  conditions,  the  old  solutions  failed  to  predict  the  time  when 
motion  is  stopped  or  the  final  shape  of  the  target. 

The  need  for  comparing  results  with  these  prime  physical  ob¬ 
servables  led  to  to  undertake  a  new  investigation  of  this  problem.  At 
this  time  specimens  of  deformed  target  plates  of  rolled  armor  steel 
became  available  frcn  Watertown  Arsenal  Laboratories ,  etched  by  Ober- 
hoffer’s  reagent.  These  yielded  values  of  deformation  under  the  pro¬ 
jected  impact  by  direct  measurement.  Even  though  relatively  few  such 
specimens  were  available,  they  constituted  an  opportunity  for  an  im¬ 
proved  theoretical  approach  to  the  problsn.  Following  suggestions  in 
the  literature  which  assumes  the  predominant  forces  in  the  deformation 
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of  a  plate  under  impact  to  consist  of  frictional  forces,  the  new  featur 
is  to  suppose  there  is  a  threshold  stress  level  or  impact  yield  const  an 
at  which  flow  starts,  and  below  which  the  t  ^  either  elastic  or 

more  simply,  rigid.  This  behavior  prevents  the  material  from  flowing 
indefinitely. 

The  analysis  of  this  non-linear  problem  is  described  in  Techni¬ 
cal  Report  no.  3>  presentation  no.  1  and  publication  no.  3-  Values  for 
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impact  yield  constants  of  210  x  lO''  psi  were  found  for  mild  steel  and 
deformations  agreed  very  well  with  those  shown  on  the  etcned  specimens. 

3*  Pe- etration  in  Fiberglass  Reinforced  Plastics  (FRP) 

This  investigation  was  undertaken  to  determine  the  basic  know¬ 
ledge  of  impact  behavior  of  a  class  of  materials  (FRP)  and  was  especial 
ly  spurred  on  by  the  needs  of  the  military  for  light  personal  armor. 

A  series  of  70  penetration  tests  were  completed  of  small  caliber  pro¬ 
jectiles  through  laminated  fiberglass  plates.  The  material,  in  proper 
combinations,  was  found  to  give  a  substantial  saving  cf  weight  ( 1U>  to 
30£)  over  steel  for  the  same  stopping  power.  A  firing  range  was  built 
and  projectile  velocities,  both  incident  ani  resiiual,  were  measured  at 
6  stations  by  aluminum  foils,  in  conjunction  with  a  Polaroid  camera  ani 
an  externally- triggered  oscilloscope .  (Publication  no.  6) 
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5*  Other  Completed  Work 

Our  work  on  viscoelastic  waves  and  further  studies  in  penetra¬ 
tion  dynamics ,  which  have  not  previously  been  reported,  are  covered  in 
detail  further  below  in  this  report . 

6.  Recent  Work  (last  6  months)  Not  Completed 

One  of  the  really  large  problems  in  dynamics  is  the  analysis 
of  stress  wave  propagation  in  multi-dimensional  bodies.  The  analytical 
solutions  to  date  apply  only  to  very  simple  geometries.  Our  success 
with  the  use  of  discrete  approaches  for  the  various  plane  and  spherical 
geometries  has  led  us  to  attempt  to  apply  this  method  to  the  afore¬ 
mentioned  problems.  To  date,  some  test  cases  hs.ve  been  solved. 

The  same  general  approach  is  also  applicable  to  a  class  of 
statical  plasticity  i  oblems.  To  date  we  have  validated  this  method 
for  some  known  problems  in  the  literature. 

The  analysis  of  flexural  travelling  waves  has  produced  solutions 
for  step  and  ramp  moment  inpvts,  but  step-shear  input  is  still  unsolved. 


Viscoelastic  Waves  with  Reflection 
for  Longitudinal  Impact  (Ch.  I-IV) 
by 

M.  L.  Wenner 


Amor  Penetration  Dynamics  (CH.  V-VIII ) 

by 


R.  Minnich 
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NOMENCLATURE 


VISCOELASTIC  WAVES 
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Mass  density  of  bar  material 
Longitudinal  coordinate 
Time 

Velocity  in  x-direction 

Cross-sectional  area  of  bar 

Tensile  stress  in  x-direction 

Tensile  strain  in  x-direction 

Pressure  at  origin 

Glassy  (fastest)  wave  speed 

Time  duration  of  stress  input 

Viscoelastic  relaxation  modulus 

Spring  constant-standard  linear  solid 

Spring  constant-standard  linear  solid 

Dashpot  viscosity  coefficient-standard  linear 
solid 

Number  of  bar  elements 
Number  of  time  elements 
Change  in  x-coordinate 
Change  in  time 

Glassy  state  relaxation  modulus 
Change  in  strain 
Length  of  bar 
Change  in  velocity 
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Longitudinal  Coordinate 
Time 

Compressive  stress  in  X-direction 

Compressive  strain  in  X-direction 

Particle  velocity 

Mass  density  of  bar  material 

Static  yield  stress  (compressive) 

Cross-sectional  area  of  bar 

Material  constant 

Overstress  exponent 

Mass  of  striking  body 

Slip  factor 

Length  of  bar 

Initial  velocity  of  striker 

Dimensionless  coordinate  (■  X/L) 

2 

Dimensionless  time(»  aQT/pDL  ) 

Dimensionless  stress  (■  o/oq ) 

Dimensionless  velocity(«  V/DL) 

Dimensionless  impact  velocity  (•  Vq/DL) 

2  2 

Dimensionless  strain  (■  oQ  e/pD  L  ) 

Dimensionless  ..'ass  factor(»  G/pAL) 

1  2 

The  quantity  k  v 

L  o 

Number  of  bar  elements 
Number  of  time  elements 


MV 


u.  ***_  ■ 


Time  at  which  unloading  begius 
Time  at  which  all  motion  stops 
Change  in  dimensionless  velocity 
Change  in  dimensionless  strain 
Change  in  dimensionless  stress 
Change  in  dimensionless  time 
"is  replaced  by" 
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Chapter  I 


INTRODUCTION 


1 . 1  On^-Dimenslonal  Impact 

Since  the  problem  of  one-dimensional  elastic  stress  vcsvcs 
in  a  slender  bar  was  solved  by  Newton  in  1685,  it  has  been  recog¬ 
nized  that  the  one-dimensional  problem  is  quite  valuable  for  de¬ 
veloping  and  checking  further  analyses  on  more  complicated  cases. 

As  our  store  of  knowledge  increases,  however,  we  find  that  the 
material  representations  which  are  most  realistic  are  usually  those 
which  also  present  the  most  difficult  mathematical  obstacles.  Thus, 
even  the  simplest  cases  of  dynamic  problems  are  difficult  to  treat 
analytically. 

It  Is  proposed  to  analyze  in  this  thesis  two  problems  of 
longitudinal  impact  of  slender  bars.  The  ultimate  concern  here  is 
to  obtain  solutions  for  realistic  materials,  since  only  in  this  way 
may  it  be  expected  that  a  contribution  of  a  quantitative  or  of  an 
engineering  nature  may  be  made. 

The  problems  to  be  analyzed  are  stress  waves  in  a  viscoelas¬ 
tic  material  and  plastic  Impact  of  a  viscoplastic  material.  The 
material  representation  for  the  viscoelastic  material  is  one  which 
can  be  directly  obtained  from  experiments,  namely,  the  relaxation 
modulus  in  tension.  The  conventional  method  cf  describing  a  material 
by  a  spring-dashpot  model  will  be  discarded.  The  law  used  herein  to 
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describe  the  viscoplastic  material  will  be  a  nonlinear  stress-strain 
iate  law. 

In  the  following,  two  problems  will  be  explained  concurrent¬ 
ly  since  many  of  their  aspects  are  quite  similar.  When  there  is  a 
dissimilarity,  separate  treatments  will  be  used. 

1 . 2  Introduction--Viscoelastic  Waves 

With  the  increasing  use  of  viscoelastic  materials  in  such 
diverse  applications  as  insulations,  binders  for  solid  fuel  rocket 
propellants,  and  structures,  there  has  occurred  an  increase  in  the 
number  of  theoretical  investigations  into  the  subject.  These  have 
been  hampered,  particularly  in  dynamic  analyses,  by  formidable 
mathematical  difficulties,  which  are  normally  overcome  by  representing 
the  material  by  means  of  spring-dashpot  models.  These  models  simplify 
the  analytical  problems  considerably,  but  at  the  expense  of  inade¬ 
quately  describing  the  materials.  In  fact,  the  commonly  used  models 
consisting  of  two  or  three  elements  exhibi'  nearly  all  of  the  change 
in  a  particular  viscoelastic  function  in  a  single  logarithmic  decade 
of  time,  whereas  experiments  show  that  actual  viscoelastic  materials 
require  at  least  seven  to  fourteen  logarithmic  decades  of  time  to 
describe  their  full  range. 

The  problem  of  one-dimensional  wave  propagation  in  visco¬ 
elastic  bars  is  among  the  simplest  of  dynamic  problems  to  state,  but 
no  analytic  solution  has  yet  been  obtained  which  does  not  depend 
in  one  way  or  another  upon  a  material  representation  of  springs  and 
dashpots.  A  general  solution  to  this  problem  would  determine  stress 
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(or  strain)  in  a  hat  as  a  function  of  time  and  position,  using 
only  the  experimentally  obtained  material  data  and  the  boundary 
conditions  on  each  end  of  the  bar.  It  is  one  of  the  purposes  of 
this  report  to  develop  an  analysis  which  yields  such  a  solution. 

In  analyzing  the  problem  of  one-dimensional  wave  propagation 
in  viscoelastic  media  the  following  assumptions  are  made: 

a)  The  bar  is  assumed  to  be  slender  enough  that 
we  have  plane  waves.  Hence,  lateral  effects  are  disre¬ 
garded,  and  the  stress  is  uniform  on  any  cross  section. 

b)  The  displacements  are  infinitesimal. 

c)  The  bar  consists  of  a  linear  viscoelastic 
material  and  thus  obeys  Boltzmann's  superposition  principle. 
This  principle,  which  forms  the  basis  of  all  linear  visco¬ 
elastic  analysis,  will  be  stated  in  Section  2.2  below. 

d)  The  bar  is  assumed  to  have  a  uniform  cross 

section. 

A  complete  solution  to  this  problem  must  show  the  response 
of  a  viscoelastic  bar  subjected  to  impulsive  loading  as  a  function 
of  time  and  position.  It  will  be  readily  adaptable  to  free  or 
fixed  ends,  and  will  show  the  reflection  of  the  waves  from  each  end. 
It  is  hoped  that  this  will  provide  a  direct  means  by  which  to 
compare  theoretical  results  with  experimental  results. 

1 . 3  Introduc tlon-Viscoplastlc  Impact 


In  order  to  adequately  describe  the  pheomenon  of  deformation 
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of  materials  in  which  plastic  deformation  occurs  at  a  high  rate,  it 
has  been  found  necessary  to  incorporate  a  constitutive  law  which 
takes  into  account  the  strain  rate.  A  simple  type  of  experiment 
is  the  axial  impact  of  a  prismatic  bar  by  a  finite  mass. 

One  of  the  laws  proposed  for  this  problem  is  a  power  law 
relation  between  strain  rate  and  dynamic  overstress.  The  advantage 
of  choosing  this  law  is  that  some  material  data  are  available  from 
test  results  and  the  law  has  shown  to  be  effective  in  predicting 
deformations  of  cantilever  beams  even  when  the  physical  constants 
had  to  be  crudely  estimated.  This  past  success  should  be  a  stimulus 
to  further  investigation  in  order  to  determine  if  the  law  in  question 
adequately  describes  different  problems.  If  this  is  the  case,  the 
theory  could  be  used  together  with  experiments,  in  order  to  determine 
the  physical  constants  more  accurately,  thereby  producing  solutions 
for  other  problems. 

The  following  assumptions  and  conditions  are  imposed  on  the 
analysis : 

a)  Uniaxial  stress  is  assumed,  and  no  lateral 
effects  are  admitted. 

b)  The  material  is  rigid  viscoplastic;  that  is, 
no  strain  increment  occurs  at  a  point  unless  the  static 
yield  stress  there  is  exceeded.  This  is  a  realistic 
assumption  if  the  plastic  deformation  at  s  point  is  much 
larger  than  the  elastic  deformation. 


c)  The  striking  mass  moves  parallel  to  the  axis 
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of  the  bar  and  after  impact,  it  sticks  to  the  end  of  the  bar. 

d)  The  constitutive  law  is  a  power  relation  be¬ 
tween  the  strain  rate  and  the  dynamic  overatress.  This  law 
will  be  stated  explicitly  in  section  2.3  below. 

The  law  to  be  used  in  this  analysis  is  a  more  general  one 
than  has  been  used  previously  in  the  problem  of  longitudinal  impact 
in  viscoplastic  rods.  These  results  should  enable  researchers  to 
more  accurately  evaluate  the  worth  of  this  particular  constitutive 
law. 

1 . A  Past  Work--Viscoelastic  Waves 

Most  of  the  analytic  investigations  accomplished  to  date 

have  used  Maxwell  solids,  Voigt  solids,  or  a  combination  of  the  two. 

The  Maxwell  model,  proposed  in  1890,  is  a  spring  in  series  with  a 

dashpot,  while  tie  Voigt  model  (1892)  is  a  spring  in  parallel  with 

a  dashpot.  These  models  have  been  used  not  only  for  stress  wave 

problems  but  also  for  vibrations  and  quasi-static  problems. 

* 

Hillier  (1]  has  used  the  Maxwell  model,  Voigt  model,  and 
two  models  using  three  elements  each  as  material  representations  for 
the  problem  of  longitudinal  sinusoidal  waves  in  a  bar.  This  is  of 
course  a  vibrations  solution  and  no  transient  effects  are  considered. 

The  transient  problem  vac  treated  by  Lee  end  Morrison  [2], 
also  using  simple  model  representations.  Laplace  transform 

* 

Numbers  in  brackets  refer  to  references  listed  in  Bibliography. 
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techniques  were  used  in  order  to  solve  the  boundary  value  problems. 

Kolsky  [3 J  performed  experiments  on  a  viscoelastic  material 
by  subjecting  a  rod  to  explosives  at  one  end  and  recording  the 
response  of  the  other  end.  An  analytic  solution  was  obtained  for 
the  same  problem  by  using  Fourier  analysis.  The  theoretical  re¬ 
sults  compared  favorably  with  the  experimental  results. 

A  partial  solution  to  the  problem  is  given  by  Bland  [4J. 

In  this  analysis  a  wave  front  expansion  is  used  which  yields  a  long 
time  solution.  We  desire,  however,  to  find  a  complete  solution 
whenever  possible. 

Morrison  (5]  has  given  integral  solutions  to  several  model 
repr ±sentat ions,  including  the  standard  linear  solid  which  will  he 
described  below.  Laplace  transform  methods  were  used,  with  a 
numerical  inversion.  Our  computer  solution  will  be  applied  to  this 
model  and  the  results  compared  to  Morrison's. 

The  most  comprehensive  solution  lo  this  problem  to  date 
has  been  given  by  Arenz  (6].  In  this  study  the  material  is  repre¬ 
sented  by  a  Kelvin  model  consisting  of  a  gla»?y  spt ing  plus  n 
Voigt  elements  in  series.  Arenz  uses  this  model  with  nine  Voigt 
elements  to  represent  the  real  part  of  the  complex  shear  compliance 
of  a  polyurethane  material.  Transform  techniques  are  used  with  two 
different  methods  of  inversion.  These  methods  were  formulated  for 
other  applications  by  Cchgpery  (7,8). 

Several  review  articles  have  recently  been  published  con¬ 
cerning  dynamic  phenomena  in  viscoelastic  materials.  Mor .  refer¬ 
enda  to  these  problems  may  be  found  in  tne  papers  by  Zverev  (9), 
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Lee  [ 1 0  J ,  and  Hopkins  (11). 

1 . 5  Past  Work--Viscoplastlc  Impact 

The  constitutive  law,  a  form  of  which  we  shall  use  in  this 
analysis,  was  first  proposed  by  Hohenemscr  and  Prager  (12].  It  is 
for  a  solid  of  the  Bingham  type  in  which  elastic  deformations  are 
neglected  and  in  which  tne  strain  rate  is  related  to  the  amount  by 
which  the  stress  at  a  point  exceeds  the  static  yield  stress.  Several 
experimenters  (13,  14,  15]  have  shown  that  this  type  of  law  satis¬ 
factorily  describes  the  plastic  impact  of  beams  and  other  structures. 

Lee  and  Wolf  (16]  made  an  investigation  of  longitudinal  im¬ 
pact  on  a  rigid-plastic  bar  in  which  the  material  was  considered  to 
be  linearly  strain-hardening  but  rate  independent.  7n  using  their 
solution  to  analyze  tests  made  by  Habib  (17],  Lee  and  Wolf  showed 
that  rate  dependence  effects  may  become  of  importance  if  nonuniform 
strain  distribution  resulting  from  plastic  wave  propagation  is  ig¬ 
nored  . 

The  linear  form  of  the  rigid-viscoplastic  lav.  has  formed  the 
basis  of  several  studies.  Sokolovskii  (18]  used  it  to  solve  several 
problems  of  plane  shear  waves  in  semi  -  inf inite  media,  and  Ting  and 
Symonds  (19]  used  it  to  analyze  the  problem  of  longitudinal  impact 
of  viscoplastic  bars.  The  same  two  authors  also  give  (20]  some 
approximate  methoas  for  the  nonlinear  case  and  compare  these  re¬ 
sults  with  both  the  linear  case  and  with  the  simplest  rigid-plastic, 
rate  independent  problem.  They  conclude  that  for  the  case  of  high 
impact  velocity  of  a  large  impact  mass  the  assumption  of  uniform 
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final  strain  is  reasonable.  We  shall  have  some  comments  on  this 
assumption  and  shall  show  its  limitations. 

Some  review  articles  which  will  provide  more  references 
are  Lee's  review  (10]  which  treats  elastic-plastic  problems,  and  an 
article  by  Cristescj  [21]  which  describes  European  contributions. 
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Chapter  II 

MATHEMATICAL  FORMULATION 

2 . 1  Peitinent  Mathematical  Relations 

The  standard  method  by  which  physical  problems  of  this 
sort  are  solved  is  to  write  down  the  appropriate  mathematical  re¬ 
lations,  combine  them,  and  then  solve  the  resulting  boundary  value 
problem.  The  approach  used  here  also  depends,  of  course,  on  a 
certain  combination  of  physical  laws,  geometrical  relations  and 
definitions.  In  particular,  what  we  need  are:  definition  of  strain,  a 
geometrical  relation  between  strain  and  velocity,  the  impulse- 
momentum  law,  and  a  law  relating  stress  to  some  other  parameters. 
These  form  the  basis  for  each  of  the  two  analyses  given. 

2. 2  Mathematical  Formulation- -Viscoelastic  Waves 

The  physical  conditions  governing  the  one-dimensional 
viscoelastic  wave  problem  is  shown  in  Figure  1.  The  origin  is 
located  at  the  left  end  of  a  bar  of  length  L.  The  right  end  may 
either  be  fixed  or  free.  The  left  end  is  subjected  to  an  impulsive 
stress  loading  at  time  zero;  this  stress  may  remain  acting  for  any 
desired  amount  of  time.  Thus,  the  boundary  condition  at  the  origin 
is 

CT(0,t)  -  -P0[H(t)  -  H(t  )  J 


i* 


(2.1) 
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where  H(t)  is  the  Heaviside  step  [unction  and  t  is  the  time  at 
which  the  stress  is  removed.  The  boundary  condition  at  the  right 
end  is 

v(L , t)  *  0  (2.2) 

if  the  right  end  is  fixed,  or 

a(L,t)  -  0  (2.3) 

if  the  end  is  free. 

The  equation  of  motion  governing  the  one-dimensional 

case  i3 

^a/Sx  »  p  c)^u/dt^  *  (2.4) 

and  the  definition  of  strain  is 

£  ■  du/dx  (2.5) 

The  above  is  standard  analysis,  and  may  be  found,  for  instance, 
in  Timoshenko  [27]. 

We  now  require  a  stress  strain  law  for  the  viscoelastic 
material.  Normally,  this  is  obtained  from  a  model  representation. 
For  instance,  the  stress  strain  law  for  the  standard  linear  solid 
shown  below  is 


( 1/E '  )d<?/dt  +  [10 


(1  +  E/E '  )de/dt  +  E[X€. 


(2.6) 
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The  standard  procedure  is  to  Insert  equations  (2.6)  and 
(2.7)  into  (2.4)  and  solve  the  resulting  partial  differential 
equation  subject  to  the  boundary  conditions,  usually  by  transform 
techniques . 

However,  the  general  stress  strain  law  for  a  linear  visco¬ 
elastic  material  is  given  by  the  Boltzmann  superposition  principle, 
which  is,  in  fact,  the  definition  of  a  linear  viscoelastic  material. 
This  is  an  integral  equation  and  as  given  by  several  authors  [4, 

22,  23]  may  take  either  of  two  equivalent  forms.  One  representation 
is 

t 

O(t)  -  €(t)E(t  -  0 )-[  € (t )  dE( t  -t)  dT.  (2.7) 

Jo  dr 


This  is  a  Riemann-Stiel t Jes  Integral,  and  when  we  apply  integration 
by  parts, 


o(0 


E(t  -  r) 


telzl 

dr 


Jt. 


(2.8) 
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There  are  corresponding  integral  relations  giving  the  strain  when 
the  stress  history  is  known.  In  (2.7)  and  (2.8)  the  function  E(t) 
is  the  relaxation  modulus  in  tension  and  it  is  Che  stress/strain 
ratio  as  determined  from  a  relaxation  test.  In  this  test  the 
material  is  suddenly  subjected  to  a  strain  which  is  subsequently 
held  constant,  and  the  stress  required  to  maintain  this  deformation 
is  measured  as  a  function  of  time. 

Inserting  one  of  equations  (2.7)  or  (2.8)  with  (2.5)  into 
(2.4)  gives  an  integro-dif ferontial  equation  which  has  not  yet  been 
solved  mathematically.  Instead  of  attacking  this  equation,  we 
propose  to  use  the  Boltzmann  superposition  principle  as  the  basis 
of  a  so-cailed  "computer  analysis"  to  solve  the  one-dimensional 
wave  problem  for  any  linear  viscoelastic  material. 

2 . 3  Mathematical  Formulation-Viacoplastic  Impact 

Figure  2  shows  a  rod  of  length  L,  fixed  to  a  rigid  wall  at 
the  right  end  and  subjected  to  the  uniform  impact  of  a  body  of  mass 
G  on  the  left  end.  The  origin  is  located  at  the  left  end  and  the 
impact  occurs  at  time  zero.  The  boundary  conditi>  at  the  left  end 
is  obtained  by  applying  the  impul ae -momentum  law  it  X  ■  0.  That  is, 

G  dV(0,T)/c*T  +  An(0,T)  -  0-  (2.9) 

The  boundary  condition  at  X  ■  l  is 

V(l,T)  -  0.  (2.10) 

The  impul  se -momeni  uir  equation  is 


cV-t/dA 


-p  £v/(*T 


(2.11) 
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and  the  equation  relating  the  velocity  to  the  strain  (where  they  are 
continuous)  is 


de/dT  -  -dv/dx. 


(2.12) 


The  equation  relating  the  strain  rate  to  the  dynamic  over¬ 
stress  is 


C*€ 

St 


c)e 

St  " 


o 


(2.13a) 


(2.13b) 


where  D  and  p  are  material  constants.  Several  consequences  of  this 
law  should  be  noticed.  First,  upon  being  loaded,  a  material  which 
obeys  equations  (2.13)  exhibits  no  deformation  until  the  stress  ex¬ 
ceeds  the  static  yield  stress.  Second,  when  a  particular  region  is 
unloaded  from  a  stress  which  is  higher  than  the  static  yield  stress 
to  a  stress  which  is  less  than  the  yield  stress,  that  region  will  move 
as  a  rigid  body. 

Following  Ting  and  Symonds  (19)  we  introduce  nondiroensional 
variables  in  order  to  simplify  the  analysis.  These  nondimensional 
quantities  are  recorded  in  the  list  of  symbols,  and  using  these 
terms,  the  equations  (2.11),  (2  12)  and  (2.13)  become 

dv/c't  ■  -t)s/c*x  (2.14) 


Sv/Sx 


-c*r|/c)t 


(2.13) 
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£r,/<^t  ■  (s  -  1)P 

it 

I8 

|  >  1 

(2.16a) 

c'Tj/c't  ■  o 

if 

1  8 

1  <l 

(2.16b) 

respectively.  In  the  case  of  p  ■  1 

the  combination  of 

equations 

(2.14),  (2.15)  and  (2.16a)  yields 

c^v/c)x^  -  dv/St  ■  0 

if 

1. 

I>1 

(2.17) 

£^s/c)x^  -  ds/dt  ■  0 

if 

|. 

Im  . 

(2.18) 

This  is  the  one -dimens iona  1  hea ;  condition  equation  and  analytical 
solutions  may  be  obtained  for  a  number  of  problems.  In  a  region 

whe  re  J  s  J  <  1, 


dTj/c't 

-  0 

(2.19a) 

c'v/c'x 

-  0 

(2.19b) 

1  2 
's/dx 

-  0. 

(2.19c) 

In  nonditnenalonal  notation  the  boundary  conditions  at  the  right  end 
becomes 


v(l,t)  -  0. 


(2.20) 


The  boundary  condition  at  x  •  0  which  was  given  by  equation  (2.9) 

is 


k  ii  .  +  i  .  o 


(2.21) 


when  combined  with  (2.16a)  using  p  ■  1.0.  In  nondimens iona 1 
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notation  the  initial  conditions  are 

v(0,0)  *  v  (2.22a) 

'  o 

v(x,0)  -  p(x,0)  -  0  (2.22b) 

s(x,0)  -  1  for  x  >  0.  (2.22c) 

The  solution  of  this  problem  is  given  by  Ting  and  Symonds  [19]  with 

solutions  to  three  similar  problems. 

It  is  proposed  in  this  report  to  give  a  means  by  which  to 
solve  the  viscoplastic  impact  problem  using  the  nonlinear  case  of 
equation  (2.13).  This  will  be  accomplished  by  means  of  a  computer 
analysis  using  the  above  relations  as  a  basis.  The  means  by  which 
these  relations  are  written  in  finite  form  and  the  method  by  which 
we  incorporate  them  into  a  computer  approach  is  given  in  the  follow¬ 
ing  chapter^ 
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COMPUTER  ANALYSIS 

3 . I  General  Computer  Approach 

The  purpose  of  this  section  is  to  give  a  brief  account  of 
the  method  by  which  both  the  viscoelastic  wave  problem  and  the 
viscoplastic  impact  problem  will  be  solved.  Instead  of  deriving 
complicated  equations  from  basic  principles  and  then  using  the 
methods  of  numerical  analysis  in  order  to  obtain  results,  we  shall 
write  down  the  physical  laws,  definitions,  and  assumptions  for  a 
particular  problem  and  use  them  in  finite  form.  This  is  called 
"computer  analysis".  Some  of  these  relations  are  given  in  this 
section,  and  the  method  by  which  we  combine  them  for  computational 
purposes  is  illustrated  in  the  following  two  sections. 

This  method  has  several  advantages  over  standard  numerical 
procedures.  First,  the  program  is  physically  meaningful  because 
the  ph>sical  laws  appear  explicitly  and  the  program  actually  follows 
the  phyrical  processes  as  they  occur.  This  fact  makes  this  method 
more  desirable  than  the  standard  procedure  in  which  equations  are 
derived  from  simple  laws  and  then  subjected  to  finite-difference 
analysis  which  usually  renders  them  unrecognizable.  Secondly,  the 
prgratr,  is  easily  adaptable  to  changes  because  often  only  one  state¬ 
ment  need  be  modified  with  ro  change  to  the  reet  of  the  program.  In 
the  use  of  a  computer  to  solve  finite-difference  equations  which 


have  been  derived  from  the  basic  principles,  the  entire  logic  must 
usually  be  chang'd  when  one  of  these  basic  laws  is  changed.  Third¬ 
ly,  the  resulting  program  is  more  efficient  because  finite  methods 
are  used  immediately  rather  than  as  a  means  to  solve  complicated 
integral  or  differential  equations. 

Since  we  are  using  a  finite  process  to  simulate  a  continuous 
process,  some  further  features  of  the  approach  must  be  explained. 
Unless  stated  otherwise  the  following  comments  apply  to  both  the 
viscoelastic  and  the  viscoplastic  programs. 

1.  The  bar  is  divided  into  a  finite  number  of  cells  and 
the  basic  equations  are  applied  to  each  of  these  cells;  the  time 
variable  is  periodically  changed  and  the  variables  such  as  stress, 
strain,  displacement,  and  velocity  are  recorded  for  each  cell  and 
for  each  time.  We  consider  here  only  bars  of  uniform  cross  section¬ 
al  area  A.  The  length  of  a  cell  is  dx,  which  is  related  to  the 
element  of  time  dt  for  the  viscoelastic  program. 

2.  It  has  been  found  that  the  order  in  wtiich  the  basic 
equations  are  stated  is  of  great  importance  and  the  order  given  in 
our  programs  is  tnc  only  one  which  has  yielded  a  solution.  The 
reason  for  this  fact  is  that  our  finite  analysis  must  necessaril> 

do  one  thing  at  a  time  whereas  m  the  actual  physical  process  sever  i; 
changes  may  occur  at  one  time.  The  stability  of  the  proceuure  is 
directly  dependent  on  the  .'•dor  of  th**  steps. 

3.  Another  factor  which  affects  stability  is  the  ceil  size 
and  the  .size  of  the  t  ins?  increment.  Detailed  comments  on  this 
problem  will  be  found  below  in  the  discussion  of  each  problem 
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4.  The  program  should  be  constructed  in  such  a  manner  that 
minor  changes  may  be  incorporated  in  both  the  data  and  the  procedure. 
If  possible,  a  single  program  should  be  able  to  solve  several  prob¬ 
lems.  For  instance,  a  program  might  have  a  feature  by  which  it  may 
solve  a  bar  with  a  fixed  or  free  end. 

5.  The  program  mechanism  should  be  constructed  so  that  the 
physical  process  is  as  closely  simulated  as  possible.  This  gives 
an  insight  into  the  problem  and  may  lead  to  further  understanding 
and  study  of  related  problems. 

Figure  3  shows  the  notation  by  which  the  cells  and  stresses 

are  labelled.  It  should  be  noted  that  the  stresses  are  shown  in  the 

positive  direction  for  the  viscoelastic  program,  but  in  the  negative 

direction  for  the  viscoplastic  program  in  which  compressive  stresses 

are  regarded  as  positive.  In  each  of  the  programs  the  index  i  runs 

from  1  to  i  ,  and  the  time  index  k  runs  from  1  to  k  . 

m’  m 

The  cores  of  the  two  problems  are  similar  and  may  be  out¬ 
lined  as  follows. 

1.  The  law  relating  stress  to  strain  is  stated  in  finite 
form  and  the  stress  in  cell  i  +  1  is  calculated  and  recorded.  For 
the  viscoelastic  problem  thi6  step  is  a  numerical  integration  while 
a  finite  form  of  equation  (2.16a)  serves  the  purpose  for  the  visco¬ 
plastic  program. 

2.  The  strain  increment  acting  on  cell  i  +  1  may  be 
written  in  terms  of  displacement  as 


dr 


i  + 


i  '  (ui  + 1  ■  ui)/dxi  + 1 


(3.1) 
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where  the  displacements  are  given  by 


ui  +  i  -  vi  +  i  dt 


(3.2) 


and 


Ui  *  Vi 


(3.3) 


Combining  (3.1),  (3.2)  and  (3.3)  yields 


dci  + 1  *  (vi  + 1  •  Vdt/dxi  +  i 


(3.4) 


which  is  the  form  we  shall  use. 

3.  A  cumulation  process  follows  in  which  the  strain  in 
cell  i  +  1  is  replaced  by  the  addition  of  de\  +  j  to  the  original 
strain  +  This  may  be  written  symbolically  as 

Ci  +  1  €i  +  l+d€i  +  l.  (3.5) 


The  symbol  «-  is  understood  to  mean  "is  replaced  by". 

4.  The  impulse -momentum  law  may  be  written  as 

^  F  dt  ■  change  in  momentum  (3.6) 

where  the  left  hand  side  of  (3.6)  is  the  sum  of  the  impulses.  We 
apply  this  law  to  cell  i  by  noticing  that  the  impulse  on  the  right 
end  of  the  cell  is  CJ^  +  j_  A  dt  an< *  the  impulse  on  the  left  end  is 
0  A  at,  where  A  is  the  cross  sectional  area  of  the  bar.  The 
change  in  momentum  of  the  cell  is  given  by  the  mass  times  the 
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increment  of  velocity, 

dv  ~  (p  A  dx^)  dv.  (3.7) 

Using  these  quantities  in  (3.6)  gives 

(p  A  dxj  dv  *  +  A  dt  -  ci  A  dt 

or,  in  the  form  which  we  shall  use, 

dv  35  +  l  "  ai^  dt/P  dxi*  0*8) 

Note  that  this  gives  an  increment  of  velocity  in  cell  i  as  a  result 
of  a  net  impulsive  force  acting  on  that  element. 

5.  In  order  to  calculate  the  actual  velocity  of  cell  i  we 
use  the  relation 

v.  «-  v.  +  dv.  (3.9) 

ii 

This  is  actually  an  integration  process  which  recalculates  the 
velocity  in  cell  i. 

6.  The  position  along  the  bar  is  calculated  by  another 
integration  process  which  we  write  as 

Xi  +  1  "  xi  +  dxi'  (3.10) 

These  six  relations,  taken  in  the  given  order,  form  the 
basis  of  each  of  the  two  solutions.  There  are  only  minor  revisions 
in  the  problems  which  will  be  explained  further  below.  Note  that 
two  of  the  relations  are  mechanical  or  physical  laws  (l  and  4), 


-  -»  v-  ■ 
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while  the  others  are  geometrical  definitions  or  integra  >.i  or  s . 

In  each  of  the  programs,  these  six  relations,  (1)  to  (6)  are 

repeatedly  calculated  for  values  of  the  index  tanging  from  1  to  i^. 

Physically,  i  is  the  number  of  elements  making  up  the  entire  length 
m 

of  the  bar.  This  repetition  for  the  index  i  is  then  nested  in 

another  repetition  procedure  which  increments  the  time  variable  by 

means  of  the  index  k.  This  index  runs  from  1  to  k  ,  where  k  Is 

m  m 

the  total  number  of  time  elements.  It  should  be  noted  that  all 
boundary  conditions  such  a*>  stress  at  x  ■  0,  momentum  interchange 
at  x  *  0  or  a  velocity  condition  at  x  -  L  are  included  in  one  or 
both  of  these  steps.  Any  initial  conditions  (as  on  the  stress 
field  for  the  viscoplastic  problem)  are  included  before-hand. 

The  above  is  a  general  sketch  of  the  method  to  be  employed. 
The  details  of  the  solution  for  each  of  the  programs  are  given  in 
the  following  sections. 

3. 2  Computer  Theory-Viscoelastic  Waves 

As  previously  stated,  a  spring-dashpot  model  of  two  to  four 
elements  is  not  sufficient  to  describe  a  realistic  viscoelastic 
material.  Although  the  number  of  elements  in  the  model  could 
theoretically  be  increased  so  that  the  model  would  represent  any 
given  material  vt:ry  closely,  this  procedure  is  not  recommended  be¬ 
cause  it  leads  to  a  numerical  analysis  scheme  of  probhibitive  length 
and  complexity.  Therefore,  we  shall  ignore  models  of  this  sort 
entirely  and  shall  use  the  Boltzmann  superposition  principle. 

Instead  of  writing  the  integrals  cited  previously  [equations 
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(2.7)  and  (2.8)]  we  shall  deduce  the  proper  form  of  the  law 
directly  from  the  definition  of  the  relaxation  modulus  and  the 
superposition  principle.  Recall  that  the  stress  at  any  time  in  a 
viscoelastic  material  due  to  a  constant  strain  suddenly  applied  at 
zero  cime  is  given  by 

O(t)  -  E(t)  €q  (3.11) 

where  e  is  the  constant  strain.  If  the  strain  had  been  applied 
o 

not  at  time  zero  but  at  time  t  »  t,  the  relation  wou' d 

0(t)  -  E(t  -  r)fo  (3.12) 

since  we  must  use  the  value  of  the  relaxation  modulus  whose  argument 

is  the  elapsed  time  since  the  strain  was  applied.  Now  let  us  apply 

(3.12)  to  a  particular  cell,  say  cell  i  +  1.  Then  the  increment 

of  stress  on  cell  i  +  1  at  time  t^  due  to  an  applied  strain  de^  ^  ^ 

at  time  t  .  .  will  be  given  by 
m  +  1 


d°i  +  1  “  d€i  +  l(tm  +  1*  X  E<S 


+  1** 


(3.13) 


In  order  to  calculate  the  stress  on  the  cell  due  to  a  series  of 
successively  applied  strains  the  Boltzmann  superposition  principle 
is  used.  This  states  that  we  may  use  (3.13)  to  compute  the  quantity 
dai  +  1  for  **c^  8traln  applied  and  add  these  quantities  in  order 
to  obtain  the  value  of  the  total  stress.  Thus,  the  total  stress 
is  given  by 

°t  ♦  1  ■  x  E"k  •  V  +  dCl  *  l(t2>  *  E(tk  *  C2> 

+  .  •  •  +  de.  +  x(tk)  x  E(t  -  0)  (3.14) 


tends  to 


UO. 


If  we  take  the  limit  as  k  tends  to  infinity  and  df^  +  ^ 
zero  we  obtain  the  familiar  representation  of  equation  (2.8): 

o(t)  »  f  E(t  -  t)  de (t)  dT.  (2.8) 

Jo  dT 


The  equation  (3.14)  is  not  quite  the  form  which  is  used  in  the 

analysis.  This  is  because  the  propagation  of  stress  waves  is  a 

continuous  process  and  (3.14)  was  written  for  a  finite  number  of 

constant  valued  strain  increments.  Therefore,  in  order  to  obtain 

the  correct  value  of  wave-front  stress  and  to  more  accurately 

describe  the  continuous  process,  the  terms  E(t,  -  t  ,  ,)  are  re- 

placed  by  Eft.  -  -r(t  +  t  ,  ,)].  That  is,  the  strain  is  regarded 
r  '  1  k  2  m  m  +  1  ' 

as  being  applied  at  a  time  halfway  between  the  times  t  and  t  ,  . . 

It  is  imperative  to  do  this  in  order  that  a  correct  value  of  the 

wave  front  stress  Is  obtained.  If  E(t.  -  t  )  were  used,  the  wave 

km  ’ 

front  stress  would  not  be  relaxed  at  all  and  if  E(t,  -  t  ,  ,)  were 

k  m  +  1 

used,  the  stress  would  be  relaxed  too  greatly. 

A  different  form  of  the  integral  equation  can  be  obtained 
by  applying  an  integration  by  parts  to  equation  (2.8).  The  result 
is 

a(t)  -  € ( t)  E(t  -  0)  -  f  e(t)  dE( t  -  r)  dr.  (2.7) 

Jo  dt 

This  representation  ia  attractive  because  it  quite  strikingly  shows 
two  physical  phenomena.  The  first  term  on  the  right  represents  the 
stress  acting  on  the  element  as  a  result  of  the  strain  which  has 
deformed  it  at  the  instant  considered.  This  part  is  completely 


analogous  to  the  case  of  elastic  waves  in  which  the  function  E(t) 
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would  be  a  constant.  The  integral,  however,  represents  in  our 
problem  the  process  of  relaxation  in  which  the  stress  which  has 
been  acting  on  the  element  for  some  finite  time  is  reduced  in  value 
because  of  the  nature  of  the  relaxation  modulus.  Ip  order  to  apply 
the  analysis  of  fhis  report  to  the  integral,  we  again  write  in¬ 
crements  of  stress  and  add  these  increments.  In  this  case  the 
quantity  d(j^  +  is  calculated  by 


d0i  +  1 


€i  +  l(tm  +  1J fE(tk  +  1  '  C 


m  +  1 


)  -  E(t 


k  4  1 


-  t 


m 


•t  2 


)] 


(3.15) 


Inherent  in  this  equation  is  the  assumption  that  the  strain  on 

element  i  +  1  is  constant  from  time  t  to  time  t  ,  .  This 

m  m  +  i 

equation  is  calculated  for  the  entire  strain  history  of  the  element 
and  the  stress  is  thus  integrated.  The  complete  the  calculation 
of  equation  (2.7),  the  "elastic"  part  of  the  stress  is  calculated 
by  means  of 


doi  +  i  -  +  i<'k  +  i>  *<*  -  °>-  <3-w 

It  was  found  that  an  effective  way  to  increase  stability  in 
this  program  was  to  compute  the  integral  in  two  places  in  the 
procedure  of  the  problem.  Therefore,  the  order  of  the  laws  given 
in  section  3.1  was  modified  as  follows.  (1)  The  first  calculation 
is  the  evaluation  of  the  integral  as  described  above.  (2)  The 
increment  of  strain  de .  ,  (t,)  is  calcuiated  and  added  to  the 

1  +  1  K 


strain  on  cell  i  +  1  at  the  last  time  index,  k  -  1.  Thus  we  cal¬ 
culate  the  strain  €,  (t  ),  which  had  not  as  yet  been  available. 
(3)  This  newly  calculated  strain  is  then  used  to  obtain  the  elastic 
part  of  the  stress  by  means  of  equation  (3.16).  (4)  The  Integral 

is  again  evaluated,  but  now  the  strain  €.  ^  ,(t  ,  . )  has  been  cal- 

culated  for  m  +  1  ■  k,  whereas  it  had  not  been  available  for  the 
first  integral  calculation.  (5)  The  stress  ai  +  ^  is  now  calcu¬ 
lated  by  adding  the  elastic  part  of  the  stress  to  the  average  of 
the  two  integral  calculations. 

The  reason  that  it  was  found  necessary  to  incorporate  this 
somewhat  artificial  device  is  that  if  it  were  not  used,  the  integral 
would  be  evaluated  either  by  using  all  of  the  strain  +  j ( t^) ,  or 
by  using  none  of  it.  Either  of  these  two  alternatives  may  lead  to 
serious  errors.  For  instance,  a  stress  may  be  calculated  which  is 
higher  in  magnitude  than  the  input  stress,  or  a  change  of  signs  may 
occur.  Either  of  these  errors  tend  to  be  magnified  as  the  calcu¬ 
lation  proceeds,  rendering  the  results  meaningless.  In  addition, 
our  method  is  probably  closer  than  either  alternative  to  the  actual 
physical  process  in  which  several  changes  occur  simultaneously. 

The  only  other  exceptional  feature  of  this  program  is  that 
the  elements  of  time  need  not  be  equal.  Since  elements  of  time  are 
related  to  elements  of  distance  (cell  size)  by  the  relation 

dxi  "  cg  x  dti  (3.17) 

it  follows  that  the  cell  sizes  may  likewise  be  unequal.  If  these 
sizes  are  unequal,  an  interpolation  procedure  is  required  in  order 


that  the  proper  values  of  the  relaxation  modulus  are  available.  A 
linear  interpolation  is  used  which  introduces  some  errors  into  the 
glassy  and  transition  regions  of  the  modulus;  these  are  not  serious 
if  the  time  elements  are  small  enough. 

Two  schemes  have  thus  been  devised  for  the  calculation  of 
viscoelastic  stress  waves.  The  remainder  of  what  follows  in  this 
section  applies  to  both  representations. 

In  order  to  allow  these  schemes  to  be  used  for  bars  with 
fixed  or  free  ends,  an  "end  factor"  i c  employed.  This  is  simply  a 
number,  either  one  or  zero,  which  is  multiplied  with  the  stress 
°im  +  i*  ^  the  bar  is  *ree  aC  tbe  r*8bc  end,  EF  -  0.0,  and  if  the 
right  end  is  fixed,  EF  =  1.0.  Thus,  the  stress  at  the  end  of  the 
bar  is  made  to  be  zero  for  a  free  end,  and  is  undislurbed  for  a 
fixed  end.  These  conditions  lead,  by  means  of  the  impulse-momentum 
law,  to  the  proper  type  of  reflection  at  the  right  end. 

The  left  end  condition  is  decided  by  means  of  a  test  in 

the  tune  repetition  loop,  but  outside  the  position  loop.  This  test 

applies  a  given  constant  stress,  either  tensile  or  compressive,  to 

the  first  cell  if  the  time  is  less  than  t  and  applies  no  stress  if 

P 

the  time  is  greater  than  t  .  It  would  be  a  simple  matter  to  apply 

P 

any  given  stress  at  any  time  to  the  bar,  but  since  this  would  add 
essentially  no  new  information  or  new  understanding  to  the  problem, 
it  was  not  done  in  our  program. 

The  basic  laws  as  they  appear  in  the  double  integral  program 
are  shown  ir.  figjre  h.  Note  that  this  diagram  is  merely  the  logical 
"skeleton"  of  the  program. 


FIG.  4  VISCOELASTIC  COMPUTER  SCHEME  USING  TWO 
BOLTZMANN  SUPERPOSITIONS 


3  -  }  Computer  Theory- Vi  sc oplast  ic  Impact 


This  program  was  construct*,  u  along  the  general  lines  given 
in  section  3.1.  A  minor  difference  is  that  the  strain  increment  is 
calculated  and  integrated  before  the  stress-strain  rate  law  is 
applied.  This  order  was  chosen  because  the  stress-strain  rat.*  law 
is  not  used  when  there  is  no  deformation  occurring,  whereas  the 
strain  is  always  calculated.  Thus,  the  order  i  merely  a  matter  of 
convenience,  since  this  particular  choice  does  not  alter  the  logic. 
The  basis  of  the  scheme  is  outline  in  Figure  5. 

In  the  computer  procedure  we  write  the  stress-strain  rate 
law,  equation  (2.13)  as 


st  + 1  ■  1-°  f  i(vi  -  \  +  i>/dx'l/p  if  + 1 


(3.13a) 


si  +  1  "  1'° 


if  V  <v.  A  , 
1  l  +  l 


(3.18b) 


This  representation  uses  a  somewhat  different  criterion  than  does 
equation  (2.13).  In  (3.18)  it  is  the  velocity  that  determines 
whether  or  not  deformation  occurs  whereas  in  (2.13)  the  stress  is 
the  determining  factor.  This  difference  arises  because  the  initial 
input  to  the  bar  is  a  velocity  on  the  left  end.  This,  of  course, 
gives  rise  to  stresses,  but  we  follow  the  diffusion  of  velocity 
through  the  bar  and  use  it  as  a  "deformation  criterion' . 

The  ms  by  •jhich  this  deformation  decision  is  made  is 
the  use  of  "slip  factors".  This  is  a  factor  (called  si.  m  the 
program)  which  takes  the  value  of  1.0  *f  cells  :  and  i  -  1  ire 


deforming  relative:  to  each  other,  and  the  value  0.0  if  the  cells 
i  and  1  -  1  are  moving  with  the  same  velocity.  We  thus  have 


si.  =  1.0 

l 

for 

Vi  -  1 

4  \ 

si.  =  0.0 

for 

v. 

=  V. 

1  1  -  1  1 


The  dynamics  of  cell  number  i  is  thus  determined  by  the 

relative  values  of  the  slip  factors  si.  and  si,  ,  .  We  compare 

l  l  +  l 

these  two  slip  factors  by  computing  their  difference,  si .  +  ^  -  si.. 
The  results  of  this  subtraction  may  lead  to  any  one  of  three 
possibilities;  that  is,  the  difference  may  be  zero,  positive,  or 
negative.  There  are  actually,  then,  four  dynamic  conditions: 

(Al.  si.  =  si.  .  «  0.0) 

A.  si  -  si.  -  0  1 

1  (A2.  si.  «  si,  ,  -  1.0) 

i  i  +  1 

B.  si.  -  si.  >  0  (si,  =  1.0,  si.  *  0) 

i  +  T  l  i  +  1  '  l 

C.  si,  .  -  si.  <  0  (si.  =  0.0,  s'  »  1.0) 

This  method  was  devised  by  Minnich  and  Davids  [25]  for 
another  application,  and  the  details  of  their  method  are  similar  in 
some  respects  to  those  used  here.  However,  the  method  will  be 
given  explicitly  here  since  the  physical  conditions  of  the  two 
applications  are  different  in  a  number  of  circumstances. 
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Condition  Al  si .  =  si .  .  ,  =  0,0 

- -  i  l+l 


In  this  case,  since  each  slip  factor  is  zero,  the  three 
velocities  shown  above  are  equal  and  there  is  no  net  force  on  cell 
i.  Therefore,  its  velocity  increment  is  zero  and,  if  we  write  the 
impulse -momentum  law  in  the  form 

dv  «=  (s^  x  8l.  *  s,  +  j  x  slj  +  pdt/dx  (3.19) 

we  may  use  (3.19)  to  calculate  dv. 

Condition  A2  si ,  ■  si .  .  ,  ■  1.0 

.  i  l  +  l 


k9> 


Now  deformation  occurs  at  both  ends  of  the  i  -  th  cell  and 
therefore  there  are  stresses  on  each  end  of  it.  Thus,  equation 
(3.19)  may  again  be  used. 

Condition  B  si.  «  0.0.  si.  .  =  1.0 

- -  L  x  +  i 


Physically,  this  means  that  the  material  immediately  to  the 
left  of  the  i  -  th  cell  is  not  deforming,  while  the  material  to  the 
right  is  deforming.  Thus  the  region  in  question  is  in  the  "un¬ 
loading"  process  in  which  the  boundary  between  rigid  and  deforming 
material  is  moving  to  the  right.  Ting  and  Symonds[19]  prove  that 
this  unloading  must  start  at  the  impact  end  of  the  bar,  and  since 
this  is  a  diffusion  phenomenon,  once  a  region  has  unloaded,  it  will 
not  deform  again.  We  thus  use  a  form  of  the  impulse -momentum  law 
in  which  the  mass  of  the  striker  must  be  caken  into  account.  This 
is 

dv  =  -s.  dt/k.  (3.20) 

i+l 


This  velocity  increment  is  added  to  the  first  i  cells. 
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Condition  C  si .  *  1.0,  si ,  ,  =  0.0 

-  i  l  +  1 


In  this  instance  the  material  to  the  left  of  the  i  -  th 
cell  is  deforming  while  the  material  to  the  right  is  rigid.  This 
corresponds  to  a  loading  process  where  the  deformation  field  is 
moving  to  the  right.  Since  s^  +  ^  must  be  1.0  and  s^  must  be 
greater  than  1.0,  deformation  will  be  initiated  to  the  right  of 
cell  i.  Therefore,  8^  +  ^  is  set  equal  to  one  and  equation  (3.19) 
is  used. 

The  structure  of  the  logic  for  this  program  is  indicated 
in  the  schematic  diagram,  Figure  6. 


i 


FIG.  6  LOGIC  OF  SLIP  FACTORS 
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Chapter  IV 

RESULTS  OF  CALCULATION  AND  DISCUSSION 

4. 1  Standard  Linear  Solid  Material 

Two  computer  programs  based  on  the  theory  given  in  sections 
2.2,  3.1  and  3.2  were  written  in  the  Fortran  language.  These 
programs  are  given  in  their  entirety  in  Appendices  A  and  B. 

In  order  to  check  the  validity  of  the  programs,  runs  were 
made  using  the  standard  linear  solid  to  represent  the  material.  The 
values  of  the  system  parameters  were  chosen  as  follows: 

E  -  E'  -  1.0 

1/H  -  1.0 

p  -1.0 

p  -  -1 .0 
ro 

The  relaxation  modulus  in  this  case  is 

E(t)  -  1  +  e-t  (4.1) 

Actually  hese  values  were  chosen  in  order  to  facilitate  a  comparison 
with  both  Morrison  [5]  and  Arenz  [6].  As  Figure  7  shows,  the  in¬ 
tegral  solution  by  Morrison,  the  Laplace  transform  method  of  Arenz, 
and  the  computer  solutions  agree  very  well. 


COMPARISON  OF  SOLUTIONS  FOR  STANDA 
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As  a  point  of  interest,  we  show  in  Figure  8  the  response  of 
a  bar  of  the  same  material  and  at  the  same  location.  The  differ¬ 
ence  this  time  is  that  the  bar  has  a  finite  length  (L  ■  2.828),  md 
the  right  end  is  fixed.  A  constant  unit  stress  is  applied  at  the 
origin.  The  step  discontinuities  in  the  response  are  due  to  the 
reflections  from  the  ends  of  the  bar. 

Several  facts  should  be  noted  from  the  diagrams.  First, 
the  response  of  the  bar  takes  place  in  less  than  two  decades  of 
log  time,  thus  bearing  out  what  was  stated  previously:  that  the 
standard  linear  solid  is  a  fictitious  material.  Also,  if  the 
response  of  two  bar  locations  are  plotted,  the  slope  becomes  less 
steep,  as  we  should  expect. 

4. 2  Realistic  Viscoelastic  Material 

We  now  apply  the  program  to  a  realistic  material.  Visco¬ 
elastic  data  was  taken  from  a  thesis  by  Arenz  [6]  for  a  polyurethane 
synthetic  rubber,  a  low  modulus  polymer.  The  relaxation  modulus  is 
shown  in  Figure  9.  There  are  some  grap..s  given  in  Arenz 's  work 
showing  stress  wave  behavior  as  calculated  by  an  approximate  Laplace 
inversion  technique.  The  single  integral  program  was  applied  to 
this  material  and  a  comparison  of  results  is  shown  in  Figure  10. 

The  high  frequency  response  as  calculated  by  Arenz  is  slower  than 
our  data.  This  seems  to  be  true  generally.  Also,  Arenz  ootained 
some  oscillation  in  the  high  frequency  response,  supposedly  due  to 
alternate  reinforcement  and  interference  of  waves  of  differing 
frequency  and  therefore  differing  speeds  of  propagation.  No  such 
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occurrence  was  obtained  in  our  analysis,  and  indeed  it  seems  likely 
that  the  oscillation  reported  by  Aren::  is  not  a  genuine  physical 
fact,  but  a  result  of  some  mathematical  approximation.  The 
transition  region  of  the  response  is  quite  well  matched  for  the 
two  solutions  but  some  divergence  is  apparent  at  large  time.  It 
is  believed  that  this  is  due  to  an  inaccurate  low  frequency  material 
representation  in  this  analysis. 

The  double  integral  program  was  applied  to  another,  similar 
material  in  reference  [24],  This  program  does  not  seem  to  operate 
as  effectively  as  the  single  integral  program,  and  there  was  some 
scattering  of  results.  Figure  11  shows  the  response  of  this 
material  (Hysol  8705)  at  two  positions  along  the  bar  x  •  1.77  inches-  ^ 
and  x  «  3.29  inches. 

It  was  found  that  several  factors  could  give  rise  to  in¬ 
stability  in  either  of  the  programs.  First,  the  time  element  must 
be  chosen  small  enough  so  that  enormous  changes  in  the  relaxation 
modulus  do  not  take  place.  This  factor  is  far  more  critical  for 

the  double  integral  program  than  for  the  single  integral  solution. 

-3 

The  time  interval  was  taken  as  3.0  x  10  sec.  for  the  results  in 

Figure  10  and  as  1.0  x  10  ^  sec.  for  those  in  Figure  11.  Secondly, 

since  each  of  the  programs  incorporates  a  linear  interpolation,  the 

material  data  must  be  given  to  the  program  at  quite  a  few  points. 

Normally,  the  data  was  given  at  log  time  increments  of  0.1  and  even 

this  is  not  enough  for  response  calculations  at  positions  where  the 

*6 

glassy  wave  speed  arrival  time  is  less  than  10  sec. 


The  double  integral  program  was  so  written  that  variable 


I 
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elements  ol  time  could  be  chosen.  It  was  found,  however,  thit 
actually  using  variable  time  elements  leads  to  serious  errors  in 
most  cases.  Generally,  it  may  be  said  that  the  time  increments 
should  be  decrease'4  in  size  as  the  time  increases.  This,  however, 
is  not  an  advantage  over  taking  equal  time  increments,  insofar  as 
required  computer  time  or  representation  of  material  are  concerned. 

It  may  be  said  that  for  all  cases  tested,  the  single  integral 
program  outperformed  the  double  integral  program  in  every  way. 
Furthermore,  it  is  more  efficient  than  the  double  integral  program. 

4. 3  Viscoplastic  Impact 

A  program  incorporating  the  theory  of  sections  2.3,  3.1  and 

3.3  was  written  in  Fortran.  This  program  is  shown  in  Appendix  C. 

The  case  of  the  overstress  exponent  equal  to  unity  was 
performed  first,  and  the  results  were  compared  with  those  of  Ting 
and  Symonds  [19].  Figure  12  shows  a  comparison  of  our  stress  cal¬ 
culations  with  those  of  Ting  and  Symonds.  This  plot  is  the  quantity 
(8  -  1)/vq  versus  the  dimensionless  distance  x.  The  calculation  is 

for  values  of  k  ■  1.0  and  v  ■  1.0.  Calculations  for  other  values 

o 

of  k  and  vq  showed  similar  agreement  with  Ting  and  Symonds.  Figure 

13  shows  a  plot  of  dimensionless  strain,  t),  divided  by  r)^  (which  is 

defined  as  the  uniform  strain  which  could  absorb  the  initial  kinetic 

energy  at  stress  s  ■  1)  versus  the  dimensionless  distance.  The 

2 

quantity  is  equal  to  one-half  the  product  k  vq  .  Again  agreement 

was  as  close  for  other  values  of  k  and  v  .  In  Figures  12  and  13  the 

o 

values  of  t  and  t^  are  the  times  at  which  the  striker  stops  moving 
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FIG.  12  COMPARISON  OF  STRESS  RESPONSE  OF  VISCO¬ 
PLASTIC  IMPACT  -  OVERSTRESS  EXP.  OF  2.0 


VISCOPLASTIC  IMPACT 


and  the  velocity  in  the  entire  bar  vanishes,  respectively. 

The  case  of  overstress  exponent  equal  to  unity  is  quite  un¬ 
realistic.  In  engineering  situations,  we  would  require  an  analysis 
using  the  correct  value  of  p.  For  mild  steels,  p  *  5,  and  for 
aluminum  alloys,  p  =  4;  therefore,  the  analysis  for  p  ■  1  does  not 
give  a  quantitative  answer  to  the  impact  problem. 

Accordingly,  we  introduced  the  value  of  p  «*  4  into  the 
program.  Thus  the  stress-strain  rate  law  becomes 

Si  +  1  “  1  +  (vi  ‘  vi  +  l/dx)Q  (4>2) 

where  Q  =  1/p.  No  serious  difficulties  were  encountered  as  long 
as  the  value  of  the  time  increment  was  kept  sufficiently  small.  For 
instance,  for  dx  =*  0.050,  the  program  operated  satisfactorily  for 
dt  =  0.001,  but  became  unstable  for  dt  =  0.01.  This  arises  because 
the  computer  performs  one  operation  at  a  time.  If  the  time  increment 
is  too  large  relative  to  the  distance  increment,  one  parameter  may 
accumulate  errors.  The  velocity  field,  for  instance,  may  reverse 
directions,  or  the  strains  become  impossibly  high. 

The  results  of  the  calculations  for  p  ■  4.0  are  shown  in 
Figures  14  and  15,  where  they  are  compared  to  the  solutions  for 
p  *  1.0.  As  could  be  seen  from  the  general  law,  the  strains  are 
larger  for  this  case.  Also,  the  stresses  are  lower  initially  but 
increase  to  a  higher  value  at  times  t  *  0.1  and  t  -  0.36.  The  time 
for  complete  cessation  of  motion  for  the  case  p  «  1.0  was  t  ■  0.770. 
At  this  time  there  was  still  plastic  deformation  occurring  in  the 
oar  when  p  ■  4.0.  It  was  also  noticed  that  the  velocity  of  the 
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X  (DIMENSIONLESS  DISTANCE) 

FIG.  14  COMPARISON  OF  VISCOPLASTIC  STRESS  SOLUTIONS 
FOR  p  »  I.O,  p  *  4.0 
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FIG.  15  COMPARISON  OF  VISCOPLASTIC  STRAIN  SOLUTIONS 
FOR  p*-  I.O,  p  *  4.0 
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striking  mass  slows  down  much  taster  for  p  ^  1  thin  for  p  -  4.  For 

instance,  with  vq  *  1.0,  k  -  1.0  the  striker  velocity  at  t  ■*  0.360 

was  0.31  for  p  ■  1.0  and  0  50  for  p  ■  4;  at  t  ■  0.600,  the  velocity 

was  0.07  for  p  •  1.0  and  0.25  for  p  ■  4. 


SUMMARY  ANJ  CONCLUSIONS 


1. .  1  Summary 

The  problems  of  longitudinal  impact  are  important  bacause 
of  their  relative  simplicity.  Their  study  indicates  the  pertinent 
physical  laws  of  a  problem  and  often  indicates  the  direction  in 
which  further  research  should  be  directed.  Even  more  important, 
they  provide  a  simple  and  direct  method  by  which  a  particular  law 
may  be  experimentally  verified.  The  two  proolems  presented  here 
belong  to  that  class  of  problems  which  have  received  considerable 
attention  in  recent  years.  They  are:  viscoelastic  waves  in  a 
longitudinal  bar  and  viscoplastic  impact  of  a  longitudinal  bar. 

Viscoelastic  investigations  typically  begin  with  a  model 
representation  of  the  material.  Son::,  of  these  studies  are  described 
in  the  first  Chapter.  When  simple  models  are  used,  however,  the 
material  is  highly  fictitious.  If  a  model  is  used  which  does  repre 
sent  a  viscoelastic  material,  the  solution  usually  involves  an  ex¬ 
tremely  difficult  numerical  progrrm.  In  this  work  a  finite  numeric* 
scheme  has  been  devised,  using  the  Boltzmann  superposition  principle 
as  toe  s.ress  strain  law.  Spring-dashpot  models  have  beeft  eliminate 
altogether,  and  the  actual  material  data  is  used  in  graphical  form. 
This  method  has  been  shown  to  solve  problems  represented  by  woopis 
as  well  as  problems  represented  by  more  realistic  materials. 
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The  current  situation  for  plastic  impact  of  bars  may  not  be 
so  easily  summarized.  There  exist  many  mathematical  models  for  these 
phenomena,  with  varying  degrees  of  complexity.  We  have  chosen  to 
analyze  the  case  of  plastic  impact  of  a  finite  mass  on  a  bar  of 
length  L,  where  the  material  is  rigid-viscoplastic.  That  is,  the 
bar  does  not  deform  until  the  stress  at  a  point  exceeds  the  static 
yield  strength.  When  it  does  deform,  it  does  so  according  to  a  power 
law  relation  between  the  rate  of  strain  and  the  amount  by  which  the 
stress  exceeds  the  static  yield  stress.  The  solution  has  been  shown 
to  agree  very  well  with  existing  analytical  results  using  the  linear 
law.  The  nonlinear  case  has  also  been  solved,  and  a  great  difference 
is  shown  between  the  linear  and  nonlinear  cases.  In  addition,  the 
final  strain  for  the  nonlinear  case  has  been  shown  to  differ  greatly 
from  that  obtained  by  assuming  uniform  strain,  when  the  impact  mass 
and  velocity  are  small.  Our  solution  also  yields  the  value  of 
stress,  strain  and  velocity  at  any  point  of  the  bar  at  any  time,  in¬ 
stead  of  just  the  final  strain. 

A«  2  Conclusions 

The  model  representation  of  viscoelastic  materials  is  in¬ 
adequate  to  describe  the  phenomenon  of  stress  waves.  The  definition 
of  a  linear  viscoelastic  material  is  the  Boltzmann  superposition 
principle  and  this  should  be  used  to  calculate  any  short  time  effects. 
The  response  of  a  realistic  viscoelastic  material  takes  place  ever 
a  large  number  of  decades  of  log  time.  This  indicates  that  phenomena 
occurring  in  material  which  is  more  rigid  than  that  used  here  will 
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require  a  great  deal  of  time  to  reach  equilibrium.  We  conclude 
that  the  study  of  viscoelastic  problems  other  than  stress  waves 
should  also  use  the  Boltzmann  principle. 

The  linear  law  of  viscoplastic  impact  does  not  give  quanti¬ 
tatively  correct  results  for  common  materials  such  as  steels  and 
aluminum  alloys.  The  nonlinear  law  also  differs  greatly  from  the 
simplifying  assumption  of  uniform  f  iin  cases  of  low  impact  mass 
and  velocity.  This  is  an  important  example  for  experimenters,  since 
it  would  be  easier  to  conduct  a  test  with  small  parameters  than 
with  very  large  ones.  Also,  the  nonlinear  law  extends  the  time  of 
the  problem;  since  this  analysis  gives  the  complete  stress,  strain 
and  velocity  distributions  in  the  bar  at  any  time,  tests  could  be 
made  to  check  all  these  quantities  for  any  period  of  time. 

A . 3  Suggestions  for  Further  Study 

The  method  of  viscoelastic  wave  analysis  presented  here 
should  be  used  in  an  attempt  to  solve  other  problems  of  more  direct 
engineering  and  research  value.  The  problems  of  two  dimensional 
waves  and  of  long  time  duration,  complicated  geometry  and  with 
accompanying  creep  are  examples  of  other  engineering  applications. 

In  the  area  of  research,  a  program  of  this  nature  might  be  used  in 
reverse  to  calculate  material  data  with  given  stress  wave  response. 
The  interesting  Fourier  analysis  of  Kolsky  [3]  in  which  he  calcu¬ 
lated  the  stress  wave  response  to  an  explosive  discharge  at  one  end 
of  a  bar  is  a  potential  check  on  our  method.  This  study  is  particu¬ 
larly  valuable  because  experimental  data  is  also  given.  Finally, 
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rhe  problem  of  impact  by  a  finite  mass  would  be  a  valuable  extension 
of  this  problem. 

The  most  immediate  and  important  use  of  the  viscoplastic 
program  would  be  to  compare  it  with  extensive  experimental  data 
If  justified,  it  could  be  immediately  applied  to  other  geometries 
and  structures.  If  the  law  is  found  lacking  in  any  way,  it  might 
be  combined  with  strain  hardening  effects  in  order  to  decide  what 
type  of  constitutive  equations  are  most  applicable  for  certain 
problems.  This  could  then  be  used  as  an  aid  in  designing  and  in¬ 
terpreting  experiments.  Eventually,  criterions  for  failure  by 
various  means  could  be  added. 
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APPENDIX  A 


COMPILE  RON  Fortran 

C  VISCOELASTIC  WAVES  IN  BAR 

BOLTZMANN  SUPERPOSITION  PRINCIPLE 
DOUBLE  INTEGRAL  PROGRAM 
LIST  OF  SYMBOLS 

X  -  LONGITUDINAL  coordinate 

T  «  TIME 

S(I)  «  TENSILE  STRESS  IN  X-DIRECTlON 
E ( I » K )  »  TENSILE  STRAIN  IN  X-DIRECTION 

V!X,T)  ■  VELOCITY  IN  X-OIRECTIQN 
AREA  «  CROSS-SECTIONAL  AREA  OF  BAR 
RHO  ■  MASS  DENSITY  OF  BAR  MATERIAL 
CG  ■  GLASSY  ( FASTEST  )  WAVE  SPEED 

TP  «  TIME  DURATION  OF  STRESS  INPUT 

PO  -  PRESSURE  AT  ORIGIN 

EMU)  «  VISCOELASTIC  RELAXATION  MODULUS 

T1(M)  -  INTERMEDIATE  TIME  ( «T ( K* 1 ) -T ( M*l> ) 

EMO  •  GLASSY  STATE  RELAXATION  MODULUS 

EMl(M)  •  FM  EVALUATED  AT  TIME'TI (M) 

TL  •  LOG  (TIME) 

IM  «  NUMBER  OF  BAR  ELEMENTS 

KM  «  NUMBER  OF  TIME  ELEMENTS 

OX  -  CHANGE  IN  X-COORDINATE 

DT  -  CHANGt  IN  TIME 

DE(I.K)  ■  CHANGE  IN  STRAIN 

DV  •  CHANGE  IN  VELOCITY 

DIMENSION  E< 22.1201 . DEI  22, 1201  • T I  200 » .DTI  200) .0X1200) .XI 200) 
DIMENSION  S ( 200 ) . V ( 200 ) . I  DENT ( 16 ) »EM! 200) » TL I  200 ) 

DIMENSION  T 1 1  200 ). EMI  1 200) 

1  READ  801.  I0ENT 

2  READ  802. RHO, TP. PO 

3  READ  802. AREA. EF 
A  READ  804,1  KM, IB. IE 
5  READ  803.  EH* 1 ) 

801  FORMAT! 16A5I 

802  FORMAT  (5F10.0) 

803  FORMAT  (E10.4I 

804  FORMAT  (8110) 

EMO-EM! 1) 

CG  ■  SORTF  <  EMO/RHO) 

T ( 1 1-0*0 
00  101  N-l.KM 
READ  805»EM!NM)  »TL ! N+l )  .A 
803  FORMAT  (E10.4»F10, 4,110) 

T  (NM )  -10.0**TL(N4l  ) 

DT !  N ) <*T ! N*1 I-T!N) 

DX(N)»CG*0T(N) 

IF  (A)  20.101.20 
101  CONTINUE 
GO  TO  IT 
20  Nl-N*l 

00  104  I-N1.KM 
EM! 1*1 )»EM!N1 ) 

OT!  n-DT(N) 

ox!  n-cr*«or  <  I  ) 

T(I*l)»T!l>*nT!N) 

104  CONTINUE 
17  EMO  ■  EM! 1 ) 

PRINT  901.  I0ENT 
PRINT  902. RHO. TP, PO 
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APPENDIX  A  (conLinued) 


PRINT  903. AREA 
PRINT  904,  I M »KM 
PRINT  90S 

908  FORMAT  ( 1 H 1 *?0hTCN5lLf  MODULUS  DATA/) 

DO  105  K«1*!M 

PRINT  909  *  TL ( K ) *  EMI K ) 

909  FORMAT  (1H  .F10.2.5X.F8.0) 

105  CONTINUE 

DO  102  K* 1 » KM 

WRITE  TAPE  2*TIK),ISII  )*I»1*K) 

XU  J-0.0 

IF  IK-1)  41*41*40 
C  INTERPOLATION  PROCEDURE 

40  DO  110  N«  1  *  K 
TlINl-TIK+D-TIN+l) 

DO  111  J-l.KM 

IFU1  (N)~T(  J+l  M  10*30*111 
111  CONTINUE 

30  EM1(N)«(EM( J+l >-EM( J)  >*(  THN1-TI  JII/(TU*1  )-T  I  J>  )+EMI J) 
110  CONTINUE 

C  STEP  PRESSURE  INPUT.  LEFT  END 

41  IF  I  T I  K 1  “TP  I  11.12*12 

11  P-PO 

CiO  TO  13 

12  P-3.0 

C  PROPAGATION  PROCEDURE 

13  SU1—P/AREA 
00  103  I-1.IM 
Kl-K-1 

IF  IK-1)  27,27*25 
25  A-0.0 

C  ANELAST1C  PART.  BOLTZMANN  SUPERPOS I TION 
32  00  107  M-1.K1 

A-A  +  E ( I >1 *M+1 )»(EM1(M)-EM1IM-U  > ) 

107  CONTINUE 

C  DEFINITION  OF  STRAIN 

DEI  l  +  l  .KI-IVI 1+1) -VII  H  +  DTIK1/DXI  1+1  ) 

ei  i+i  ,K)-Eim»Ki  i  ♦oe  ( i+i  *k> 

C  STRESS-STRAIN  LAW.  ELASTIC  PART 

B»E 1 1 ♦ 1 »K ) *EM0 
C-0.0 

DO  109  M-1.K1 

C'C*EII+liM*l )*IEM1(M)-EM1 IM+1 ) > 

109  CONTINUE 

SI  1+1 ) -0.5*1  A+C)+B 
IF  (I-IMI  27,14,14 

14  $  ( I  +1 )  »$U  +1  > *EF 

C  IMPULSE-MOMENTUM  LAW 

27  OV-ISI 1  +  11 ‘AREA- SI  I ) *ARE A ) +0T I K I / ( RhG*aREA#OX( 1 1 ) 

VII ) - VI l )+DV 
XI I+l )-X(I )+DXl! ) 

103  CONTINUE 
102  CONTINUE 
REWIND  2 

PRINT  905,  (XII ) ,I'!B»!E) 

905  FORMAT  I  IMS*  nx»E10»4) 

00  106  K-l.KM 

READ  TAPE  2»T(K).*SU)*I»1*K) 

PRINT  90».7fK) ,ISI 1 >. J-ie,I£) 

906  FORMAT  I  IMS, ElO.fc.  IX . 19F6. 3 ) 
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106  CONTINUE 
500  5  TOP 

901  FORMAT!  1H1  , /.OX  .25HVI  SCOELASTIC  WAVES  IN  BAR/ 1H0  •  16' 

902  FORMAT! 


112H0MA5S  DENSITY 

• 

F10.1 *10H 

// 

232H  PULSE  DURATION 

m 

• 

F10.1.10H 

// 

332H  PULSE  INTENSITY 

■ 

• 

F10.1. 10H 

/) 

903  FORMAT! 

132H0AREA  Or  BAR 

■ 

t 

FI0.1.10M 

/) 

904  FORMAT! 

132H0NUMBER  OF  CELLS 

9 

• 

13.  // 

232H  NUMBER  OF  TIME  ELEMENTS 

• 

• 

I 3/1H1  1 

END 

! 

i 
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AITKNDIX  B 


’  !  C  r.  ■  .  .  ■  -  I.  T  r:  •_ 

r  v  1  '  ••  •  t  f  a  *»' 

C  mm:  t z".w  :  ».  rppKiPif 

r  MNGl.g  i  NT  f  ’ju  At  S'-./'m.AM 

(.1ST  Of  SYMBOLS 

X  *  longitudinal  coordinate 

T  *  i  i  ME 

SU>  '  TENSILE  STRESS  IN  X-DIRECTION 

Ell  *  K  )  »  TENSILE  STRAIN  IN  X-DIRECTION 
V(X.T)  «  VELOCITY  IN  X-OIRECTION 

AREA  -  CROSS-SECTIONAL  AREA  OF  BAR 

RHQ  ■  MASS  DENSITY  OF  BAR  MATERIAL 
CG  »  GLASSY  (FASTEST!  WAVE  SPEED 

T 0  •  TIME  DURATION  OF  STRESS  INPUT 

PO  »  PRESSURE  AT  ORIGIN 

EM(X)  -  VISCOELASTIC  RELAXATION  MODULUS 

EMO  ■  GL A$$w  STATE  RELAXATION  MODULUS 

T1IN)  -  INTERMEDIATE  TIME 

EMI  ( M )  -  EM  EVALUATEO  AT  TIME-T1IM) 

TL  ■  LOO  (TIME) 

! M  *  NUMBER  OF  BAR  ELEMENTS 

XM  «  NUMBER  OF  TIME  ELEMENTS 

DX  -  CHANGE  IN  X-COORDINATE 

DT  ■  CH,  NGC  IN  TIME 

DEI  It XI  -  CHANGE  IN  STRAIN 

DV  -  CHANGE  IN  velocity 

DIMENSION  El  22. 1201  .DE ( 22. 120) » T < 200 ) »TTI 20G ) 
DIMENSION  SI2  00) .VC200I .IDENTI16I .EMI  200 » »TLC 200) 
DIMENSION  T1 12001 .EMI (200) 

1  READ  801.  IDENT 

2  READ  802.RH0.TP.P0 

3  READ  802. AREA. EF 
A  READ  80A.IM.KM.IB.IE 
5  READ  803.EMI 1»,0T 

801  FORMAT! 16A5I 

802  FORMAT  (5F10.0) 

803  FORMAT  (2E10.AI 
80A  FORMAT  (8110) 

EM0»EM( 1 ) 

CG  «  SORTF  (EMO/RHO) 

OX«CG*OT 

T(n»o.o 

17  EMO  •  EMC  1  ) 

TTin-0.0 
00  101  J-1.200 
READ  809.TL!JM).EM(JM),A 
809  FORMAT  (F10.A.E10«A»I10I 
IF  (A)  90.90.91 

90  TTt  JM)"10.0*»TL(JM) 

T( JaI )»TI J)*OT 

101  CONTINUE 

91  KMM-J-I 
PRINT  901.10ENT 
PRINT  902.RH0.TP.P0 
PRINT  903. AREA. OY 
PRINT  90A.  IH.KM 
PRINT  908 

908  FORMAT  (1H1 .20HTENSILE  MODULUS  OATA/J 
00  10S  X-l.lM 
PRINT  909. TT IK) .EMIK ) 
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..LVKN.'H X  i'.  (continued) 


coo  (  )  M  iH0.4.SX,F|n.?l 

los  ccntinuf 

no  102  K-l.KM 

W R 1  IE  TAPE  2  •  T ( K ) • ( SI  1  )  • 1 - 1  •  K ) 

X-0.0 

C  » NTEPPOl AT  ION  PROCEDURE 

< 1-K-l 

IF  U-  1>  41.41  * 4 0 

40  T1U  >«T IK )-0.5aDT 
IF  IK-2)  42  *  42  *  TO 

TO  00  110  N*2*K1 
T 1 1 N)«  T) |N-i 1-OT 

110  CONTINUE 

42  DO  111  N-l.K 
00  112  J- 1 »  KMM 
IF  I  THNI-TTl  JM>»  30.30*  112 
112  CONTINUE 

30  EMI  IN)  •  (EMI  J.  1  i -EM  (  J)  1*1  T 1 1 N » -  T  T  C  J>  )/ITTI  JM  I-TTUI  »*£MI  J) 

111  CONTINUE 

C  STEP  PRESSURE  INPUT.  LEFT  END 

41  IF  (Tirj-tF)  11,12,12 

11  P«PO 

60  TO  13 

12  P-0.0 

C  PROPAGATION  PROCEDURE 

13  SI  1 )«-P/APEA 
DO  103  I-l.IM 
K 1-K-l 

IF  U-l)  27.27,50 
C  DEFINITION  OF  STRAIN 

50  DEI  I  ♦  1  .<  >  -  « V I  hll-Vl  m«DT/DX 
EII*l.TI*EII*l.iUl*OE(l*l.KI 
C  ANFLASTIC  PART.  BOLTZMANN  SUPERPOSITION 

CN-0.0 

DO  10<J  1  •  K  1 

CN-CN-DEI  m,M*ll»lEMl  |M|) 

109  CONTINUE 
S 11*1 »  * Cn 
f F  I  I  -  1  M  I  27,14,14 

14  SI !M ) -SI !♦! >*EF 

C  IMPjlSE -MOMENTUM  LAW 

27  OV-  ISIIM)»«EA*SH  »  •  AREA  I  »DT  /  I  RMO«*AR£A*DX  » 

VI 11 -VI  I »*DV 
X-X-OX 

103  CONTINUE 
102  CONTINUE 
REWIND  2 

do  ios  r»:.*w 

READ  TAPE  ?,T  IX  I  ,  I  SID  .  I  *■  I  ,X  ) 

PRINT  906,TU).(sm*l«iB,lE) 

VOo  FORMAT  I 1MS.E10.4. 1X.19F6.3I 
lOfc  CONTINUE 
530  STOP 

•01  FORMAT! 1M1 .40X.25MV1SCOELASTIC  WAVES  IN  BAR/ 1H0, 14AS » 

02  FORMAT! 


132MOMASS  DENS1 Ty 

•  .  F10.1.10M 

// 

23?m  PULSE  DURATION 

-  .  F10.1.10M 

n 

33?M  PULSE  ‘NTFNSITY 

•  »  FI o.l.l OH 

n 

903  FORMAT! 

132MOAREA  OF  BAR 

•  .  F10.1.10M 

// 

Al’PKND  I X  B  (continued) 


?‘'2H  ?!  ZF  Of  T  !M£  ELEMENT 
'’(K  f  OPM A  T ( 

n?HONUMPEP  OF  CELLS 
?3?H  NUM6EP  of  time  elements 
L  NO 


•  E10.<**10H 

»  13.  // 

»  I3/1H1I 


n  r  n  n  nriAnAAario^A^nnAnn 
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AIT:'. NO  [ 


t\ 


c 


COMPILE  RUN  f'OPTRAN 
C  VlSCQPI.  ASTIC  IMPACT  OF  (TAPS 

F!M!Tc  MARS  IMPACT  ON  LEFT  F NO 
POWER  STRESS  -  STRAIN  RATE  LAW 

dimensionless  form 

LIST  OF  SYMBOLS 
X  «  DIMENSIONLESS  coordinate 
TT  =  DIMENSIONLESS  TIME 
S  «  D  I  Mf‘is  I  ONLESS  STRESS 
V  *  DIMENSIONLESS  VELOCITY 
VO  =  dimensionless  IMPACT  VELOCITY 
E  =  DIMENSIONLESS  STRAIN 
P  <  OVERSTRFSS  EXPONENT 
D  «  MATERIAL  CONSTANT 
XK  »  DIMENSIONLESS  MASS  FACTOR 
ETA  *  THE  QUANTITY  0. 5*X K* ( VO ) **2 *  0 
EE  «  E/ETA 
SS  -  (S-WO)VO 
SL  ■  SLIP  FACTOR 
IM  -  NUMBER  of  bar  elements 
<M  -  number  of  TIME  ELEMENTS 
DV  *  CHANGE  in  dimensionless  VELOCITY 
Dt  *  change  in  dimensionless  strain 
dt  ■  changf  in  dimensionless  time 

DS  «  CHANGE  IN  DlMr  SIGNLESS  STRESS 

DIMENSION  S( 500 1 *V( 500) *  I  DENT 1 16) *S$(500) iE ( 500 ) »DE ( 500 ) » EE ( 500) 
DIMENSION  SLI500) 

1  READ  801.  IDENT 

2  READ  802.V0.DT .DX.XK.P 

3  READ  803.IM.KM 

801  FORMAT  (  HA5) 

802  FORMAT  (5F10.5I 

803  FORMAT  (2110) 

TT=C.0 

PRINT  901. IDENT 
PRINT  90c. VO. DT. DX.XK.P 
PRINT  903.IM.KM 
CO  101  I  *  1  *  I  M 

s  m »  i  .  o 

VII) »o.o 
101  CONTINUE 
Q-1.0/P 
VII )»V0 
DO  102  <-l »KM 

WRITE  TAPE  2.TT.(S(I).I*l.IM).S(IM+ll,IV(I)»I»l.IM)»IEm.I«l.IM) 
X-0.0 

C  IMPACT  OF  FINITE  MASS  ON  LEFT  END 
IF  C 1-1  )  17,17.18 

18  DV1 ■ <-DT/XK)»( 1.0- <V( 2 )-V< 1 ) )/DX)*»Q 
IF  IVIll+DVl)  19,17,17 

19  VU  )«0.0 
GO  TO  20 

17  V  ( 1  )  -V (  1H-DV1 

20  SLI  1  )  =*0*0 

DO  103  1*1 . 1 M 

14  DEI  I  +  1 ) ■ I V ( 1  +  1 ) “VI  I  I )*DT/DX 
e 1 1 +i ) *r 1 1 +1 ) +  DE 1 1 ♦ 1 ' 

IF  (Vin-VIItl))  29,28.28 
29  VI  I+D-VI I  ) 

C  VISCOPLARTIC  STRESS  STRAIN  RATE  LAW 
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APPi-'N'DtX  C  (continued) 


28  S  f  1*1  1  =  1  *  0*  f  (VC  1  I  —V I  l+l  )  I /DX  >  **Q 
C  TE^T  FOR  SLIPPING 

2 7  IF  !V(  I  +  1 I-VI  I  I  )  11*30*30 

30  SL( 1+1  )*0.0 
GO  TO  32 

31  SI!  I  1 »  *  1  •  0 

32  IF  ( 1-1 )  11*11*36 

11  IF  tSL(*n  33  *33.12 

12  OV-O.O 
GO  TO  35 

36  IF  (SL! I+ll-SL! I >1  37,33.34 

37  SL! 1+11*1.0 
GO  TO  3? 

34  0V»  l-1.0l»SII  *'1)*0T/XK 
IF  IVUI+OV)  40*40.41 

41  vi  n*v<  i  i+ov 

GO  TO  42 
40  VI  I )-0.0 

42  DO  104  1 1*>  1  *  I 
VII 1 )«V( I  1 

104  CONTINUE 
GO  TO  15 

C  IMPULSE  MOMENTUM  LAW 

33  DV* !S( I )*Sl( I  )-S( 1  +  1 )*$L< 1+1 ) )*DT/DX 

35  vm*v<  D+DV 

15  X=X+DX 
103  CONTINUE 
T  TeTT  +DT 
102  CONTINUE 
REWIND  2 
DO  105  K-l.KM 

READ  TAPE  2  .  TT  *  I  SU)  *  l  •  1  *  IM)  ,S  ( IM+1) ,  <  VII)  .  I  ■  1 » IM)  ,  I  E  (I)  ♦  I-  1 . 1  M ) 
imi-im+i 

DO  107  1=1, IM1 

s  s  m * ( s  u  >  - 1 .  o  i  /vo 

107  CONTINUE 

PRINT  906, TT,  <SSm»I»l*lMl> 

906  FORMAT!  1HS»F5.3*1X*21F5. 2) 

105  CONTINUE 
REWIND  2 

DO  106  KB 1  *  KM 

READ  TAPE  2  *  TT  »  (  S<  I)  *  1  - 1 , 1 M ) , 5 ( I M+ 1 ) , ( V ! I > , I ■ 1  *  I M) * ( E  (  I ) ,  I  ■  1 , 1 M ) 

PRINT  907, TT.  <  V ( I) ,I«1,JM> 

907  FORMAT  ( IMS , F5 *  3  *  3 X *2 QF5 . 2 ) 

106  CONTINUE 
REWIND  2 

ETA»0.5«XK*( V0**2.C> 

DO  108  K-l.KM 

READ  TAPE  2.TT*  (  S(I» *  I  *  1 .  IM) ,S (I M+1J , I V< I) *  I -1  * IM) , < E (I) *  I ■ 1 . 1 M > 

DO  109  I-l.IM 
EE! I  )*E( I) /ETA 
109  CONTINUE 

PRINT  908.TT.  (EE  (I  )  *  l  ■  1  *  1 M  ) 

908  FORMAT  ( 1HS.F5.3 *3X *20F5 .2  I 

108  CONTINUE 
500  STOP 

901  FORMAT  <  1 H 1  ,40X , 27HV I SCOPLAST  I C  IMPACT  OF  OARS/ 1H0 * 16A5 ) 

902  FORMAT! 

1 32H0I mPACT  VELOCITY  “  *  F10.3#10H  // 

232H  ELEMENT  OF  TIME  -  ,  F10.3.10H  // 
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APPENDIX  C  (continued) 


>32H  El  CMENT  Of  OAR  LENGTH 

*  » 

F10.3. )0H 

U 

A32H  DIMENSIONLESS  MASS  VAC  TOR 

*  • 

F10. 3. 10H 

// 

532H  OVERSTRtSS  EXPONENT 

a  * 

FI  0.3.1  OH 

n 

903  EORMATT 

132H0NUM8ER  OF  BaR  ELEMENTS 

"  ♦ 

I  3  •  H 

232H  NUMBER  OF  TIMt  ELEMENTS 

■  » 

I 3/1H1  1 

END 
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CHAPTER  V 

AN  ANALYSIS  OF  ARMOR  PENETRATION  DYNAMICS 

by  R.  Minnich 

5. 1  Introduction  and  Assumptions 

A  reasonable  set  of  assumptions  which  may  be  made  In  analyzing 
the  motion  of  a  projectile  as  It  penetrates  armor  material  are: 

I)  The  projectile  Is  assumed  to  be  a  non-deforming  body  of 
arbitrary,  but  known,  geometry  and  mass. 

II)  The  projectile's  motion  during  penetration  Is  resisted  by 
a  system  of  forces  which  depend  upon  geometry,  Initial 
velocity,  and  the  material  properties  of  the  armor.  These 
forces  are  assumed  to  be  of  two  types:  the  resistance  of 
the  material  to  penetration  due  to  Its  compressive  resist¬ 
ance  and  the  Inertial  resistance  of  the  material  as  It  Is 
displaced  by  the  projectile. 

Ill)  Frictional  effects  are  neglected  at  present  but  could  be 
added. 

5.2  Derivation  of  the  Governing  Equation 

Because  frictional  effects  are  neglected,  the  resisting  force 
components  are  assumed  to  be  acting  norma)  to  the  projectile  surface. 
Figure  12  shows  the  force  component  acting  on  the  elemental  surface 
area,  dA$. 
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The  compressive  resistance  force  Is  considered  to  be  uniformly 
distributed  over  the  surface  area  of  the  projectile  tip.  This  force 
per  unit  area,  a  property  of  the  armor  material,  will  be  denoted  as  O'. 

The  inertial  resistance  is  not  uniformly  distributed  over  the 
surface  area  of  the  projectile  but  Is  dependent  upon  the  shape  of  the 
projectile.  To  find  this  force  It  Is  assumed  that  the  change  In 
kinetic  energy  of  the  armor  material  being  displaced  Is  equal  to  the 
work  done  by  the  Inertial  force  on  an  element  of  armor  material.  This 
relation  can  be  expressed  as 

df  dx  .  1/2  v  2  dm  (5. 1) 

n  n  n 

where: 

dfn  ■  normal  force  acting  on  a  differential  area  (dAs)  of 
the  projectile  surface 

dx  =  displacement  of  the  element  of  mass  of  the  armor 
n 

material  normal  to  the  projectile  surface. 

v  »  velocity  of  the  element  of  mass  in  a  direction  normal 

to  the  surface  of  the  projectile  (equal  to  the  normal 

component  of  the  projectile  velocity). 

dm  *  mass  of  the  differential  element  of  the  armor  material 

being  displaced  (equal  to  p  dx^  dAs,  where  p  (s 

the  mass  density  of  the  armor  material). 

The  substitution  of  p  dx  dAs  for  dm  and  v  cos  &  for  v 

n  * 

yields 

df  »  1/2  o  (cos20)  v2  dAs 

n 


(5-2) 


66* 


From  these  equations  Adams  and  Tsai  (II)  derive  an  equation  of 
motion  given  by 


»•* 


! 


O’  cos  8  dAs 


■n 


p  (cos^S)  v2  dAs 


(5-3) 


where  the  component  of  each  force  per  unit  aree  In  the  direction  of 
motion  It  Integrated  over  the  frontal  surface  area  of  the  projectile. 
They  then  proceed  to  solve  this  problem  In  three  parts  for  three  con** 
dttlons  which  arise,  depending  upon  the  location  of  the  projectile 
nose  in  the  armor  material.  The  three  positions  are:  the  entrance 
phase  when  the  surface  area  of  the  projectile  Increases  with  the 
depth  of  penetration,  the  phase  where  the  nose  is  completely  imbedded, 
and  the  exit  phase  where  the  area  decreases  to  zero  as  the  nose 


emerges. 
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CHAPTER  VI 

COMPUTER  ANALYSIS  OF  PENETRATION 


6. 1  I ntroduct ion 

The  theory  developed  in  Chapter  V  is  not  limited  to  the  calcu¬ 
lation  of  residual  velocities  for  complete  penetration.  Much  more 
complicated  mathematical  procedures  would  be  required  to  solve  the 
governing  equation  to  predict  depth  of  penetration  if  complete  pene¬ 
tration  does  not  occur.  To  avoid  these  limitations,  a  computer 
analysis  was  developed  which  combines  the  derived  expressions  for  the 
forces  with  the  physical  laws  governing  the  system  and  sidesteps  the 
mathematics.  The  program  has  all  the  advantages  mentioned  in  Section 
3.1.  The  resulting  program  is  general,  and  enables  one  to  predict 
the  depth  of  penetration  or  residual  velocity  depending  upon  the 
initial  conditions. 

6.2  Development  of  the  Program 

The  projectile  was  divided  into  cells  by  cross-sections  normal 
to  its  axis.  A  typical  cell  is  shown  in  Fig.  12,  its  length  being 
dxl  and  its  surface  area  being  dAs.  The  method  employed  was  to  sum 
the  forces  acting  on  each  cell  and  use  an  impulse -momentum  relation¬ 
ship  to  calculate  the  velocity  change  of  the  projectile  for  each  time 


i  nterva 1 . 
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As  was  stated  in  Chapter  V  the  forces  are  assumed  to  consist 
of  two  kinds,  an  inertial  force  and  a  compressive  resistance  force. 

Let  f j  and  denote  the  sum  of  the  components  in  the  direction 

of  motion,  of  the  compressive  force  and  the  inertia  force  respectively. 
Then  dfj  represents  the  component  in  the  direction  of  motion  acting 
on  an  individual  cell  of  surface  area  dAs  and  is  defined  as 

df  j  =  O'  cos  0j  dAs.  (6. 1) 

where  the  subscript  i  denotes  the  i-th  cell.  Similarly,  df^  is 
given  by 

df2  =  ^  p  cos^0.v2dAs.  (6-2) 

The  sum  of  fj  and  will  then  represent  the  total  force  acting 

on  the  projectile  during  a  time  dt.  This  sum  is  given  as 

f  =  f,  +  f2  (6.3) 

The  velocity  change  dv  will  then  be  calculated  using  the 
impulse-momer.tum  law.  This  is  expressed  as 

dv  *  -  fdt/H  (6.4) 

Because  the  forces  act  opposite  the  direction  of  motion,  a  minus  sign 
is  included  in  the  above  equation.  This  velocity  increment  is  then 
added  to  the  velocity  of  tha  projectile  to  obtain  the  velocity  at  a 


given  time  or 


V 


v  +  dv 


(6.5) 


The  distance  the  projectile  travels  during  each  time  interval, 
dx,  is  obtained  by  integrating  the  velocity  over  the  time  Interval  dt. 
This  is  expressed  as 


dx  =  vdt  (6.6) 

The  total  depth  of  penetration  or  the  distance  from  the  surface  of 
the  plate  where  the  initial  contact  is  made  to  the  nose  of  the  pro* 
jectile  is  the  sum  of  all  these  incremental  dx's. 

x  x  +  dx  (6.7) 

Equations  6.1  -  6.7  are  combined  in  a  program  which  calculates 
the  force  acting  on  the  projectile  during  each  time  Interval  and  finds 
the  velocity  change  during  this  interval.  This  process  begins  when 
the  projectile  first  penetrates  the  plate  and  stops  when  the  projec¬ 
tile  velocity  becomes  zero  (for  partial  penetration)  or  when  the  pro¬ 
jectile  leaves  the  plate.  We  again  add  a  00  loop  which  directs  the 
repetitive  operations  from  i  a  I  to  i  =  i  to  obtain  f  and  use 
another  00  loop  to  repeat  this  process  for  successive  time  intervals. 
The  program  then  appears  as  below. 


vdt 

definition 

-  fdt/H 

Impulse-momentum  law 

v  ♦  dv 

x  +  dx 

output 


The  above  is  the  basic  program  except  for  goometry  calculations,  which 
will  be  discussed  below,  input  statements  and  output  statements,  and 
various  tests  to  determine  if  the  plate  has  been  completely  penetrated 
or  if  the  projectile  has  stopped. 


6.3  Gtometry 

The  various  projectile  configurations  dealt  with  and  the  dimen* 
sions  which  need  to  be  specified  to  completely  define  them  are  shown 
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in  Fig.  13-  Two  quantities  Had  to  be  determined  for  each  cell  in 
order  for  the  force  calculations  to  be  carried  out,  namely  the  surface 
area,  dASj ,  and  the  cosine  of  the  angle  between  the  normal  of  the 
cell  and  the  direction  of  motion,  cos  0j.  These  quantities  were 
determined  at  the  beginning  of  the  program  and  stored  until  they  were 
needed.  Their  calculation  for  each  projectile  configuration  follows. 

6. 3» i  The  Ogive 

Fig.  14  shows  the  ogive  with  its  detailed  dimensions.  The 
quantities  yQ,  and  y  were  required  to  be  computed  for  the  cal* 
culation  of  the  cosine  and  the  surface  area  for  each  cell.  OB  repre¬ 
sents  the  entire  surface  of  the  nose  while  BH  is  the  cylindrical 
surface.  OP  is  the  axis  of  the  projectile.  dAs.  is  the  surface 
area  of  the  i-th  cell  and  dxl  is  its  length.  The  origin  of  the  x 
and  y  axes  was  taken  as  shown.  The  needed  quantities  can  now  be 
expressed  in  terms  of  the  known  quantities  for  each  cell.  That  is 


(«2  ’  xl32) 


(6.8) 


y^  *  (A2  -  (xlj  -  xlSj)2) 


(6.9) 


y  -  y3  -  Y0 


(6.10) 


dAs.  can  therefore  be  expressed  as 


dAs.  *  2  *A  ( 1 . 0  "  yQ/y j) d*  1 


(6.11) 


and  the  cosine  as 
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FIG.  13  VARIOUS  PROJECTILE  CONFIGURATIONS 
WITH  DIMENSIONS 


9*. 


cos  8.  -  (xlj  -  xlSj)/R  (6.12) 

For  the  cylindrical  portion  of  the  projectile  the  surface  area  of 
each  cell  is  constant  and  is  given  by 

dAs.  =  2n  dxl  (6. 13) 

the  cos  0.  in  this  region  is  equal  to  zero. 

6*3.2  Hemispherical  End  Cap 

The  hemispherical  end  cap  projectile  was  useful  to  study 
because  of  the  availability  of  experimental  results. 

The  cylinder  with  hemispherical  cap  is  a  special  case  of  the 
ogive  with  R  =  r ^  and  xl^  =  r^.  Therefore  the  only  dimensions 
required  to  define  it  are  Tj  and  xl.  The  surface  area  for  elements 
of  equal  width  Is  the  same  for  a  sphere  and  is  given  by 

dAs .  =  2n  r j  dx  1  (6.14) 

This  is  also  the  expression  for  dAs.  in  the  cylindrical  portion. 

The  cosine  of  the  spherical  portion  is  defined  as 

cos  0.  =  1.0  -  (xlsj/rj)  (6.15) 

and  again  the  cosine  of  the  cylindrical  portion  is  zero. 

It  should  be  pointed  out  that  a  cylinder  with  a  hemispherical 


nose  was  programmed  rather  than  a  sphere  because  it  is  more  general. 
Since  the  forces  only  depend  upon  the  surface  area  of  the  nose,  a 
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sphere  can  be  made  a  special  case  of  the  hemispherical  cylinder  by 
setting  the  length  of  the  cylindrical  portion  equal  to  2/3  of  the 
radius  of  the  sphere.  In  this  way  the  mass  will  be  equal  to  that  of 
a  sphere. 

6.3.3  Conical  End  Cap 

The  cone  is  defined  if  r  xl,  and  xlj  are  known.  The 
distance  from  the  projectile  axis  to  the  surface  at  xlsj  is 


y,  =  r}  xlSj/xl, 


(6.16) 


The  surface  area  of  the  i-th  cell  which  is  the  surface  area  of  the 
frustum  of  a  cone  i s 


The  cosine  of  a  1 1  the  cells  is  a  constant  value  equal  to 


cos  d.  =  r  j/(r  ^  +  xl 


(6.18) 


The  surface  area  and  the  cosine  for  the  cylindrical  portion  is  the 
same  as  the  other  two  cases. 

These  three  configurations  are  related  so  that  a  separate  pro¬ 


gram  for  each  is  not  required.  The  given  dimensions  read  into  the 
program  determine  which  configuration  is  being  considered.  The  details 
of  this  are  shown  as  comments  in  the  program  itself  found  in  Appendix  0. 
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CHAPTER  VI I 

RESULTS  AND  CONCLUSIONS  OF  PENETRATION  STUDY 


7* 1  Verification  of  the  Program 

The  Initial  computations  were  concerned  with  checking  the  pro* 

gram's  results  with  those  of  Adams  and  Tsai.  The  easiest  method  of 

doing  this  was  to  see  If  a  computed  curve  of  residual  velocity  vs 

Initial  velocity  agreed  with  theirs  using  the  same  initial  conditions. 

The  initial  conditions  were 

plate  material  «  polyethylene 

plate  thickness  *  0.65  inches 

■|x  j  1a 

plate  mass  density  *  0.89  x  10  lb-sec  /In 

2 

compressive  resistance  *  25,670  lb/in 

projectile  shape  >  sphere 

The  resulting  curve  from  the  program  is  shown  in  Fig.  15  as  the  curve 
for  the  spherical  projectile.  Because  it  coincided  very  closely  with 
the  Adams  and  Tsai 's  curve  there  was  no  need  to  include  both. 

All  the  work  in  (11)  was  shown  to  agree  with  experimental 
studies.  However,  because  they  conducted  the  experiments  themselves 
it  was  decided  to  check  the  program  with  other  experimental  results. 
Gupta  and  Oavlas  (17)  performed  studies  on  the  penetration  of  fiber¬ 
glass  reinforced  plastics. 

Because  the  compressive  resistance,  <s,  is  a  property  of  the 
armor  material  it  must  first  be  found.  In  order  to  find  CT  one  must 
know  all  the  conditions  of  a  ballistics  test.  If  the  initial  velocity 
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and  the  residual  velocity  are  known,  different  values  of  cr  may  be 
tried  in  the  program  with  the  intention  of  finding  one  which  yields 
the  known  residual  velocity. 

This  process  is  not  so  much  of  a  trial  and  error  effort  as  it 
might  seem.  A  curve  of  residual  velocity  vs  cr  may  be  drawn  which 
will  enable  one  to  find  the  correct  value  of  cr  from  the  known  value 
of  vr> 

The  program  was  checked  with  the  first  three  entries  in  Gupta 
and  Davids'  Table  2.  The  third  entry  which  was  for  a  plate  0.22 
Inches  thick  impacted  by  a  0.22  caliber  projectile  with  an  Initial 
velocity  of  1270  ft/sec  was  used  in  the  calculation  of  cr,  the 
intent  being  to  determine  a  value  which  gave  a  residual  velocity  of 
750  ft/sec.  This  value  was  found  to  be  127,500  lb/ln  .  The  program 
then  correctly  predicted  the  residual  velocities  for  the  different 
plate  thicknesses  found  in  (17).  The  results  are  summarized  in  Table 
3. 


Table  3 

Comparison  with  Experimental  Results  of  (17) 


No  Thickness  Initial  Residual  Velocity 

(Inches)  Velocity  ft/sec 

_ (ft/sec)  Experimental  Program 


0.22 

1270 

750 

750 

0.09 

1270 

1090 

1085 

0.13 

1270 

1000 

995 

3 
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When  the  value  of  cr  i s  found  the  program  becomes  capable  of 
predicting  results  of  different  Initial  conditions  for  any  projectile 
configuration.  It  should  be  remembered  that  this  value  Is  only  a 
property  of  that  plate  material  for  which  It  was  found. 

Because  of  the  excellent  agreement  shown  In  Table  3  with  exper- 
(mental  results  It  was  concluded  that  the  program  Is  valid. 

7.2  Review  of  Significant  Runs 

Curves  showing  residual  velocity  vs  Initial  velocity  were 
desired  for  the  other  two  projectile  configurations  (the  ogive  and  the 
cone).  They  are  shown  In  Fig.  14  also.  The  compressive  resistance 
constant  was  obtained  from  (11)  for  epoxy  and  from  (12)  for  aluminum. 

All  of  the  curves  In  Fig.  15  approach  a  straight  line.  It  can 
be  seen  that  this  straight  line  portion  has  the  same  slope  regardless 
of  plate  material  or  projectile  configuration.  This  straight  line 
portion  has  a  slope  of  unity  which  means  that  a  change  in  the  Initial 
velocity  will  produce  the  same  change  in  the  residual  velocity. 

Figures  16,  17,  and  18  show  the  force  vs  time  plots  for  the 
three  projectile  configurations.  It  is  seen  that  they  are  all  of  the 
same  basic  shape  with  the  only  differences  being  the  entrance  and  exit 
regions.  The  rapid  rise  to  the  maximum  force  which  the  sphere  exhibits 
indicates  that  this  configuration  is  not  a  very  good  penetrator  as  com¬ 
pared  to  the  ogive.  The  cone's  ability  to  penetrate  varies  from  being 
the  best  to  the  worst  penetrator  depending  upon  the  angle  the  surface 
«rea  of  its  nose  makes  with  its  axis. 
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FIG.  16  FORCE  VS.TIME  FOR  OGIVE  PROJECTILES 
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Tables  4,  5.  and  6  contain  the  results  of  all  the  runs  made 
for  ogives,  spheres,  and  cones  respectively.  These  tables  contain 
the  runs  previously  mentioned  plus  some  other  general  runs  for  other 
plate  materials  and  initial  conditions.  These  were  made  merely  to 
show  the  applicability  of  the  program. 

It  can  be  concluded  that  for  penetration  In  which  the  projec¬ 
tile  Is  net  deformed,  this  analysis  Is  valid  both  In  predicting  depth 
of  penetration  and  residual  velocities,  and  Is  very  efficient  In  com¬ 
putational  accuracy. 


*  ar 


Table 
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2000  1885  Polyethylene  O.65  Complete  Penetration 


Summary  of  Computations  for  Hemispherical  End  Cap 
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Conclusions 

The  assumptions  and  the  analytical  solution  of  ductile  hole 
enlargement  are  presented  in  Chapter  V.  The  forces  acting  on  the 
projectile  are  assumed  to  be  of  two  kinds,  the  inertial  resistance 
of  the  plate  material  and  the  compressive  resistance.  The  expressions 
for  these  forces  are  given  and  are  then  used  to  derive  the  governing 
equation.  The  equation  for  predicting  the  residual  velocity  after 
penetration  Is  also  given. 

The  next  chapter  contains  the  development  of  the  computer  pro¬ 
gram  used  to  solve  the  ductile  hole  enlargement  problem.  The  first 
section  describes  the  incorporation  of  the  expressions  for  the  forces 
into  a  program  which  will  yield  the  velocity  change  during  a  set  time 
Interval.  The  next  section  discusses  the  geometry  of  the  three  pro¬ 
jectiles  considered  and  the  method  of  dividing  them  into  cells.  The 
expressions  for  the  surface  area  of  each  cell  and  the  cosine  of  the 
angle  between  the  normal  to  the  surface  of  the  cell  and  the  horizontal 
axis  are  also  developed. 

The  results  of  the  investigation  for  ductile  hole  enlargement 
are  contained  in  Chapter  VII.  The  first  results  presented  are  those 
needed  to  verify  the  computer  program.  Both  experimental  and  theoret¬ 
ical  ballistic  data  are  compared  with  the  output  of  the  program.  The 
comparison  in  all  cases  is  excellent.  Curves  showing  force  vs  time 
and  residual  velocity  vs  Initial  velocity  for  all  the  projectile  con¬ 
figurations  are  given.  Also  contained  in  this  chapter  are  tables 
with  the  results  of  all  the  computer  runs. 
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The  ability  of  various  projectile  configurations  to  penetrate 
a  plate  is  determined,  it  is  concluded  that  a  sphere  is  a  much  poorer 
penetrator  than  an  ogive.  The  cone  varies  in  penetration  ability  from 
the  best  penetrator  to  the  worst  depending  upon  the  angle  the  surface 
of  the  nose  makes  with  the  axis. 

8.2  Suggestions  for  Further  Study 

The  most  obvious  suggest. on  for  further  work  in  the  plug  forma* 
tion  problem  would  be  to  extend  this  investigation  to  other  materials. 
In  order  for  this  task  to  be  undertaken,  photographs  must  be  obtained 
from  which  experimental  deflection  curves  can  be  made.  All  the  condi¬ 
tions  of  the  impact  must  be  known,  i.e.,  initial  velocity  of  the  pro¬ 
jectile,  mass  of  the  projectile,  plate  thickness  and  density  of  the 
plate. 

A  criterion  for  complete  penetration  would  be  a  worthwhile 
extension.  As  was  mentioned  before  the  program  does  not  contain  this 
important  aspect.  By  inspection  of  the  final  deflection  curves,  one 
can  obtain  a  fairly  accurate  guess  as  to  what  initial  conditions  will 
cause  complete  penetration;  but  specific  results  are  needed. 

Worthwhile  studies  could  be  begun  on  other  types  of  failures 
caused  by  Impact.  A  computer  program  incorporating  the  material  laws 
associated  with  scabbing  or  dishing  wowld  be  of  interest.  Some  com¬ 
puter  work  (18)  has  been  done  for  cratering. 

Also  minor  revisions  can  be  made  in  the  present  program  to  give 
radial  strains  and  strain  rates.  A  condition  which  prescribes  the 
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initial  velocity  around  a  hole  in  the  plate  might  be  of  some  value. 

The  program  for  the  penetration  of  nonmetal  lie  materials  and 
ductile  hole  enlargement  in  metals  is  fairly  complete.  As  is  shown 
the  correlation  with  experiment  is  excellent.  Some  approaches  in 
this  area  take  into  account  a  friction  force  also.  The  program  could 
very  easily  be  extended  to  include  this  If  it  were  deemed  necessary. 

It  is  believed  that  an  extension  of  this  program  can  be  used 
to  solve  problems  of  water  entry.  For  this  problem  water  is  consid¬ 
ered  incompressible  so  that  the  compressive  force  would  be  7ero. 
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APPENDIX  A 


COMPUTER  PROGRAM  FOR  PLUG  FORMATION 


COMPILE  RUN  FORTRAN 

c  Rigid-viscous  model  for  plug  formation  in  plates 

C  UNITS  IN  IN-LB-SEC 

c  v  ■  velocity  in  z-direction 

C  Fill  ■  VISCOUS  FORCE  ON  LATERAL  AREA  OF  I-TM  CELL 

c  X.MR  ■  MASS  OF  EACH  CELL 

C  R  ■  RADIAL  COORDINATE 

C  AREA  =  INSIDE  LATERAL  AREA  OF  CELL 

C  W  ■  DISPLACEMENT  IN  Z-DIRECTION 

C  S  ■  SHEAR  STRESS  ON  R-PLANE  IN  Z-DIRECTION 

C  XMOD  ■  STABILITY  MODULUS  FOR  EACH  CELL 

C  XM  ■  CUMULATIVE  MASS  OF  FIRST  I  CELLS 

C  SI  NOND I  MENS  I ONAL  SHEARING  STRESS 

C  VV  •  NONDIMENSIONAL  VELOCITY 

C  WW  •  NONDIMENSIONAL  DISPLACEMENT 

C  TT  •  NONDIMENSIONAL  TIME 

C  DEN  ■  WEIGHT  DENSITY  OF  PLATE  MATERIAL 

C  GNU  ■  COEFFICIENT  OF  VISCOSITY  OF  PLATE  MATERIAL 

C  SY  ■  IMPACT  YIELD  CONSTANT 

c  h  •  plate  thickness 

C  VB  «  INITIAL  VELOCITY  OF  PROJECTILE 

C  RO  •  RADIUS  OF  PROJECT  ILF 

C  XMB  *  MASS  OF  PROJECTILE 

C  IM  «  NUMBER  OF  CELLS 

C  KM  »  NUMBER  OF  TIME  INTERVALS 

C  DT  *  CHANGE  IN  TIME 

C  VS  ■  VELOCITY  scale  factor 

C  WS  *  DISPLACEMENT  SCALE  FACTOR 

C  SS  *  STRESS  SCALE  FACTOR 

C  OR  *  CHANGE  IN  RADIUS 

C  FMT1  ■  VARIABLE  FORMAT  FUR  PRINT  STATEMENTS 

C  RHO  ■  MASS  DENSITY  OF  PLATE  MATERIAL 

C  XMRO  «  MASS  OF  PLATE  MATERIAL  UNDER  IMPACT 

C  VO  «  I N i T  T  AL  VELOCITY  OF  PLATE  MATERIAL  UNDER  IMPACT 

DIMENSION  VI  t  00  )  tPI  1  AOI  «XMRUOO»  »RI  10DI  »  ARE  A 1 100)  »VVf  1001  • 
1WWI 100) »W( 100) • S  t 1 00  t .si ( 100) » I  DENT ( 16  I »FMTl <141 »XMOD< 100)  • 
2  SLI100) »XM(100I 

1  READ  80’ t ( JDENT ( I  I »I» 1 ,14 ) 

001  FORMAT ( 14A5 I 

2  READ  802 •  NAME  1 .NAME  2 .DEN »GNU»$Y  »H 

802  FORMATI2A5»3F10.0.FI0.4> 

3  READ  801»VB.R0. XMB 

803  FORMATI2FJ0.3.E10.4) 

4  READ  804.IM.KM.DT.VS.WS.DR.SS.K1 

804  FORMAT(2!S.2E10.4t?r10.4.1E10,4tI5) 

IF  ( I M ) 99 .99 .5 

3  READ  801.FMT1 
PI‘3.1M5927 
PH0«DEN/384, 

XMRO*° I •RQ**2#H#RH0 
VO* I VB*XMBI /IXvBtXMRO) 

REWIND  2 

PRINT  900.1 IDENTI 11.1*1.14) 

PRINT  901 «NAM£ 1 .NAME 2 .DEN.RHO.GNU.SY 

PRINT  902  .VO.RO.H.VBtXMB 

PRINT  903.IM.KM.DT. VS. DRtWG^SS 

PRINT  905 

RID-0.00 

DO  101  I  •  1  »  1 M 

AREA  111*0.0 
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*1  !♦! )tRI| DDR 

AREA!  DIDPl’RI  I*U«2.0»H 

W(I t-0.0 

WW( I l«W< 1 )/WS 

s(  n-o.o 
sn  n«s(  i  i/ss 

FUI.0,0 

IFIRIDD -ROl  2 1 *21 *22 

21  VI  D*VO 

I  FI  1-1)15  *15*16 

15  XMI1  » - C  P I *R  f  2 ) ••2*RH0*H)»IXMB/XMRO4l«0l 
XMR ( 1 DXMID 

GO  TO  23 

16  XMI I  I « I P 1 *R ( I ♦ 1 ) •*2*RHO*H|* I XM8/XMR0+1 *0 1 
XMR ( I DXMC 1 I-XMI I “1 I 

GO  TO  23 

22  VU  1*0.0 

XMR  C  I  > ■  P I  *(R  ( IdRI  1*1 »  »»IR|  IM  I -A 1 1  1  |*RMO*M 
XM(  I  )  *XM (  1-1  DXMRl  I  I 

23  vvc  n*v<  i  j/vs 

C  MODULUS  TEST  FOR  STABILITY 

XMODI I ) *GNU*ARE  A ( !♦! DOT /(XMR  I 1)»0R» 

PRINT  906*1  *XMR< I ) *XM( 1 ) «AREA( 1 ) • XMODI I  I 
101  CONTINUE 
T  *0.0 
PRINT  904 
00  103 
TT=T/DT 

PRINT  FMTl.TT* I VV I  I  I  *  I *1 . IM) 

WRITE  TAPE  2.TT.IWWI  D*I*l.!M)»tSltlDl*l»lM) 

00  106  1*1.  IM 
W(  I  )«W(  I  DVI  n»OT 
WWI I ) =W( I J/WS 
106  CONTINUE 
SL I  D-0.0 
00  102  I  *  1  * IM 
I F  I  K  1 )  40  *40  *4 1 

40  IFU-D42.42.43 

42  PRINT  909 

43  PRINT  910.1  .VIItiSID.FlIliSK  I )  *0V 

C  R1GID-PLAST IC-Vl SCOUS  stress-strain  LAW 

41  SUM  )«-SY*GNU«!V!  Dl  1-VI  ID  /OR 
SI  (  IM  >•$<  I  +  ll/SS 

F|I*1  LSI  D1DAPEAIDD 
C  TEST  for  SLIPPING 

1FIVI  I  1-VI  011  124.24.25 

24  SLU*D«C.O 
60  TO  26 

25  SLl  IM  1*1.0 

26  IFISL  U  l-SLUM)  )2T.  28.2V 

27  IFIFID1  )-Fl  I)«XM(I)/XM(  I -D  130.31*31 

30  SLU  1-1*0 
GO  TO  2d 

31  OV*  FI  IMI/XHI  1 1 •OT 
l  F I  V I ID0VI60. 60*61 

6)  VII)*VI!)*OV 
GO  TO  62 
60  Vll 1*0.0 
62  DO  106  11*1.1 
VI ID-VI  D 
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vvi  ii  »*vi  m/vs 

108  CONTINUE 
GO  TO  35 

29  IFIFI  1  +  1MFI  IH32.32.33 

32  VU  1-0.0 
DV-0.0 
GO  TO  35 

33  5l(  1+11*1*0 

C  IMPULSE  MOMENTUM  LAW 

OV-IFI  !♦!  MSLI  I+l  I -FI  I  |  *  SL.  <  I  II/XMRI I  »«0T 
SLI 1*1 1-0.0 
GO  TO  34 

C  IMPULSE  MOMENTUM  LAW 

28  OV- (F(I*])*SL (l+l) -F(l)* SLI!) * /XMR ( l J *0T 
3*  I F  *  V « I  MDV)50.50»51 
51  VI  IMVI  I  MOV 
GO  TO  35 
50  VII 1-0.0 
35  VVI I  I -V  I  I »/VS 

102  CONTINUE 
T.T+OT 

103  CONTINUE 
PRINT  907 
REWIND  2 

DO  104  K-l.KM 

READ  TAPE  2 »TT *  I WWI I  I  *  I * 1 » 1M> ♦ I  SI  I  I ) • 1-1 » IM) 

PRINT  FMI1.TT.I WW lll.l-l.IM) 

104  CONTINUE 
REWIND  2 
PRINT  908 

DO  105  <»1  .KM 

RF  AD  TAPf  ?.TT.IWWII»»1«1»IM).IS1I||.I»1»IMI 
PRINT  FmT1.TT.IS1  I  I  I  » I • 2  *  I M I 

105  CONTINUE 
GC  TO  1 

99  STOP 

900  FORMAT t 1HI »14A5////// > 

901  FORMAT! 1M  .35X.31HVI SCO-PLASTIC  ANALYSIS  OF  PLATE  // 

132M  maT r  RIAL  ■  *2 A5 / / 

232H  wfIGMT  DENSITY  DEN  •  .F1C.5.VM  Lb/IN+*3  U 

3 si*  RhC«DEN/G  PmO  •  » E 1 0 * 4  » 1 feH  LB-SEC**2/ IN*  +  4  // 

COEFFICIENT  OF  VISCOSITY  GNU  ■  *F 10.3.13M  LB-SEC / I N**2  // 

S3?m  impact  v I F i D  CONSTANT  SY  -  »E 10.4  »9m  LB/IN+*2  I 

902  F0pMAT l 1M0.31MINI TIAl  VELOCITY  VO  -  .E10.3.7M  IN/SEC  ft 

13jh  RADIUS  OVER  lmICm  *'0  ACTS  RO  •  «Fio.3»3M  IN  // 

232H  thickness  of  plate  h  •  .fio.3.3m  in  // 

332h  VrLOCITY  oE  BULLET  VB  •  .EI0.3.7H  1N/SEC  // 

432M  “ASS  OF  PullET  *MB  •  »£10.3.13M  LB-SEC++2/IN  ) 

9G3  FORMAT  I 1M0.31MTME  NUMBER  OF  CELLS  IM  •  .110  // 

132H  Nu"-‘tf<  'JE  timl  intervals  AM  •  .110  it 

232M  C mANGI  IN  TIML  DT  •  .E10.4.4H  SEC  // 

332«  VELOCITY  SCALE  FACTOR  VS  •  .110. 4, 7M  IN/SIC  // 

432M  CHANGE  IN  RADIUS  DR  •  .F10.3.3M  IN  // 

5*?M  i)lS*»:  ACCMlsT  SCALE  FACTOR  WS  •  .FJ0.4.3H  IN  // 

63?H  STRESS  SCALf  FACTOR  SS  •  .E10.4.9M  LB/IN**2  1 

904  FORMAT l 1M1 »4hT |Ml ,5PX.?0MVrL0C I T»  PROPAGATION  //) 

903  FORMAT! 1H1 .ahCELL .8*. 9HCELL  MAES.  DX.IOmIOTAL  *ASS.9X .4MAREA* 12* t 
1  7MM0UOLUSI 

908  FORmaTIIH  » I  3 .8* »f 10.4 .8* *E I  0.4 .4* »F 1 0.4.8X.F10.4 ) 

907  FORMAT! IMl.fcMTlME .50*. 12H0ISPLACLMENT  //> 
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908  fOP.MAT  I  1H1  .4HT  IHC  #55X*6HSTRLi,S  //) 

909  FORMAT  I  1H  *1X,1HI  »)L2X#AHV(  II  *16X#4HSI  l»*16X*4HFUi  ?9X# 

1  2HSL  *15X  »2H0V ) 

910  FORMAT ( 1H  •  1 2 *5X ,F 15.4 ,5X ,F 15.4 *5X .F 15.4 »5X,F3. 1 *5X ,F 15,4  I 
END 


APFENO  I X  B 


COMPUTER  PROGRAM  FOR  ARMOR  PENETRATION 


multifile  PUN 
compile  run  fortran 

C  MMNPRCGRAM 

C  PROJECTILE  ARBITRARY  SHAPE.  NON-DEFORMING 

C  COMPRESSIVE  AND  INERTIAL  RESISTANCE  BY  TARGET 

C  SUBROUTINES  BULLET  AND  VELO  REQUIRED 

C  IF  SPHERE  SET  P’Rl  AND  XL3=Rl 

C  IF  CONE  SET  R=1000.0  AND  XL3=XL1 

C  SYMBOL  TABLE  (UNITS  IN  IN-lB-SEC  SYSTEM  I 

C  BRNO  ■  MASS  DENSITY  OF  PROJECTILE 

C  R  ■  RADIUS  OF  CIRCULAR  ARC 

C  Rl  -  RADIUS  OF  CYLINDRICAL  PORTION  OF  PROJECTILE 

C  XL  TOTAL  LENGTH  OF  PROJECTILE 

C  XLO  *  HOP.  D I  ST  ANC  t  DFT.  9  AND  TIP  OF  PROJECTILE 

C  XL  1  =  LENGTH  OF  OROJFrTlLE  NOSE 

C  X L 2  *  LENGTH  OF  CYLINDRICAL  PART  OF  PROJECT ILS 

C  XL 3  *  LENGTH  of  the  StMlCGRD 

C  DXL  *  WIDTH  OF  A  CELL 

C  H  •  PL  ML  THICKNESS 

C  TRHO  «  MAES  DENSITY  OF  TARGET  MATERIAL 

C  SIG  -  COMPRESSIVE  STRENGTH  OF  TARGET  MATERIAL 

C  DT  «  TIME  INTERVAL 

C  FS  FORCE  SCALE  FACTOR 

C  VS  VELOCITY  SCALE  FACTOR 

C  TS  -  TIME  SCALE  FACTOR 

C  VI  •=  INITIAL  VELOCITY  OF  PROJECTILE 

C  XLS  =  DISTANCE  )F  A  CELL  CRCM  TIP  OF  PROJECTILE 

C  DAS  =  SURFACE  mREA  OF  FACH  CELL 

C  ASD  *  SURFACE  AREA  OF  CELLS  IN  CYLINDRICAL  PART  OF  PROJ* 

C  COS  *  CUS.  OF  ANG.  BF T *  HUR,  DIA.  AND  NORMAL  ON  SURFACE 

C  SIN  =>  SOP T  F  (1.0-C05**2> 

C  COT  =  CCS/SIN 

C  DM  =  PROJECTILE  MASS 

C  IM  =  NCP'UER  OF  CELLS 

C  T  «  T I  ME 

C  V  r  VFIOCITY  OF  PROJECTILE 

C  X  =  DISTANCE  DFT.  TIP  OF  PROJ.  AND  FRONT  EDGE  OF  PLATE 

C  DX  DISTANCE  PROJECTILE  TRAVELS  IN  TIME  DT 

c  XP  *  DrPTH  OF  PTNETRAT ION 

c  VR  RESIDUAL  VELOCITY  OF  PROJECTILE 

c  fi  «  total  force  due  to  compressive  resistance  of  target 

c  F2  total  FORCE  dul  to  inertial  resistance  of  target 

c  F  =  SUM  OF  FI  AN"1  F  2 

c  YO  »  DISTANCE  BET.  PROJ.  AXIS  AND  HOR.  DIAMETER 

C  Y  «  RADIUS  OF  ANY  SECTION  OF  PROJECTILE 

C  Y3  «  VERTICAL  DISTANCE  OF  SURFACE  FROM  HOR.  DIAMETER 

C  VOL  I  =  VOLUME  OF  NOSE  OF  PROJECTILE 

C  VOL2  *  VOLUME  OF  CYLINDRICAL  PORTION  OF  PROJECTILE 

C  VOL  ■  VOLUME  OF  PROJECTILE  EQUAL  TO  VOL1+VOL2 

r  tt  •  mondimensional  time 

r  VV  «  NONDIMENSIONAL  VELOCITY 

C  FF  NONDIMENSIONAL  FORCE 

1  DIMENSION  XLSI200).  DAS(?00>.  COS(200).  SlNi2Q0)»  COTI200I* 

1 IDENT ( 16) »  FMTI ( 16  I 

2  Rc AD  800. ( IDENT ( I ! »  I  *  1.16  ) 

800  FORMAT  (16  AS ) 

3  REAL  806.  BNAMtl .0NAKE2 »BRHQ 
805  FORMAT  (  2A5.E10.2  I 

A  R  F  AD  8 1 0  i  R.R1 .XL3.XL.DXL 
R.n  FORMAT ( 5F 10.0) 
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5  READ  015*  TNAME1 »TNAME2*H.TRH0.S1G 
015  FORMAT C2A5.3E10.2 ) 

6  READ  020 *  DT»FS*VS»IN*NRUN 
020  FORMAT! 3E10.2.21 10) 

7  READ  03O»FMT 1 
830  FORMAT! 16A5 I 

00  400  Nc 1 .NRUN 

8  READ  025. VI 
825  FORMAT ( E 1 0*2 ) 

I  FIN- 1 ) 30.30.40 

30  PRINT  900.  (IDENT  <I).I  *  1*16  ) 

900  FORMAT ( lHl *40x .33HRESI0UAL  VELOCITY  IN  A  PROJECTILE  / 

1  1H0.25X.16A5////) 

CALL  O'JLLET  (R.R1  .XL . XL  1 »XL2 »XL3  »DXL»IN*IM»PAS  *COS*PM  tBRHOI 

SCALF  FACTORS 
TK  r  [)T 

PRINT  905  .BNAME1 »0NAMC2 . T NAME  1 ♦ TNAME2 »BRHO»OM*TRHO»S1G*H 


905  FORMAT  ( 26H0MAT!;RIAL  OF  PROJECTILE  « 

1  26 H  MATERIAL  OF  TARGET  * 

2  26 H  MASS  DEN. OF  8-MAT. BRHO  ■ 

3  26H  MASS  OF  PROJECTILE  BM  = 

4  26H  MASS  DEN. OF  T-MAT.TRHO  = 

5  26 H  COMP. RES  1 5 .OF  T-MAT.S1G- 

6  26H  THICKNESS  OF  TARGET  H  = 

PRINT  910. Xl.  .XI  1  .XL2.XL3  .R.R1  .DXL.IM 

910  FORMAT  (75HOTOTAL  LENGTH  OF  CELL  XL 

1  26H  LENGTH  Oh  VAR. SECT.  XL1 

2  26H  LENGTH  OF  CON. SECT.  XL2 

3  26H  LENGTH  OF  SEMiCHORD  XL3 

it  2  6H  RADIUS  OF  CIRC.  ARC  R 


.1X.2A5  // 

. 1X.2A5  // 

♦  E13.4.16H  LB-.,EC*«2/1N«*4// 
•E13.4.13H  LB-SEC«*2/IN// 
.E13.4.16H  LB-SEC**2/IN**4// 
.E13.4.9H  LB/IN**2  // 

. F 1 0  #  5  » 7H  INCHES  /) 

«  .F10.5.7H  INCHES  // 
e  , F10.5 »7H  INCHES  // 


. F10» 5  » 7H 
.F10.5.7H 


INCHES  // 
INCHES  // 


26H  RADIUS  OF  CON. SECT. 

26H  length  of  the  cell 

26H  NUMBER  OF  CELLS 


.F10.5.7H  INCHES  // 
.F10.5.7H  INCHES  // 
.L13.4.7H  INCHES  // 


7  26H  NUMBER  OF  CELLS  IM  «  »  110  /) 

PRINT  950.  VS.FS.DT 

950  FORMA T I 26H0VEL .  SCALE  FACTOR  VS  -  .E13.4.7H  IN/SEC  // 

1  ?6H  FORCE  SCALE  FACTOR  FS  =  .E13.4.4H  LB  // 

2  26 H  T I mf  INCRFMFNT  TS  =  OT  *  .E13.4.4H  SEC  /) 

40  PRINT  915. VI 

915  FORMATI26H1 INITIAL  VELOCITY  VI  =  .E13.4.7H  FT/SEC  /) 

C.LL  VELO  ( XL »OXL .COS » DAS *BM. IM.V I . DT »H .5 1 G* TRHO.FS * VS * 

1  TS.  1  T .X.XP.VV.FF.XLl »R1 .KM) 

REWIND  2 
Print  920 

920  FORMAT! 1H1.4HTIML.9X,11HPLNLTRAT10N.9X.8HTRAVLL-X*9X. 

18HVELOCIT Y.9X.5HFORCE  //// ) 

DO  400  K  =  1 .KM 

READ  TAPE  2.TT  .XP.X.VV.FF 

PRINT  FMT1 .TT .XP.X.VV.FF 

400  CONTINUF 

1000  STOP 
END 


COMPILE  RON  FOPTRAN 

5UBROUT INE  BULLET (R.R1 .XL  »XL1 . XL 2 » XL3 »DXL *  IN » IM»DAS »CO$*BM *BRHO » 
DIMENSION  XL SC  200 ) *DAS < 2C0) .COS! 200 ) .  S I N < 200 >  »COT ! 200) 9 V ( 200 ) 

PI  =  3.14159 
XLS! 1 )  “  0.0 
AST.  *  2  » 0  *  P  I  <  R  1 »  D  X  L 

TEST  TO  DETERMINE  IF  PPOJECTILE  IS  A  CONE 
IFIR-1000. 0)20.30. 30 
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C  LOGIC  FOR  OGIVE  AND  HEMISPHERICAL  NOSE 
20  YO  «  SQRTF  ( R** 2-XL 3**2 ) 

XLO  «  SQRTF ( R»*2- I Rl+YO) **2 ) 

XL1  *  XL 3“XLO 
XL2  *  XL-XL1 
VOL1  ■  0.0 
V0L2  ■  PI*R1»*2»XL2 
DO  101  I  «  1 • I N 
IF  (XLS( I J-XL1 )  10.10.11 

C  NON  UNIFORM  PORTION  OF  PROJECTILE 

10  Y 3  =  SQRTF(R*»2-(XL3-XLS( 111**2) 

DAS ( 1 )  =  2.0*Pi»R*( I.0-Y0/Y3)*DXL 
COSH)  •  «XL3-XLS( I »)/R 

SIN(l)  »  Y3/R 

COT  <  I  )  «  COS(  I  )  ,'S IN(  I  ) 

Y(  I)  *  Y3-YO 

0V0L1 *  PI *DXL*  < Y( 1 ) »*2+Y ( I ) *COT< I )*DXL+<COTI I »  )*»2*DXL**2/3#0) 

VOL  1  *  VOL l+DVOLl 
GO  TO  14 

11  IF  (XLS(i)-XL)  13.13.12 

12  IM  =  1-1 
GO  TO  501 

C  UNIFORM  PORTION 

13  D AS ( I )  =  ASD 
COS ( I )  =  0.0 

14  XL$ ( I ♦ 1 )  =  XLS(I)+DXL 

101  CONTINUE 

501  VOL  *  VOL 1+VOL2 
BM  =  BRHO*VOL 
GO  TO  400 

C  LOGIC  FOR  CONICAL  NOSE 

30  Y ( 1  )  =  0.0 
XL  1  =  XL3 
XL2  =  XL  -  XL  1 
DO  102  I  =  l.IN 
XLS1  1  +  1  )  =  XLS  (  I  )+DXl. 

IF  ( XL  S (  I  I  -  X  L 1  >  40.40.41 

40  Y (  1+1  I  =R1#XLS< l  +  l ) /XL  1 

HAS  C  I  )  *  PI  *<Y!I)  4  YU  +  l))*SQRTF(DXL»*2+IYm-Y(I  +  in**2) 
COS(I)  *  R 1 /SQRTF (R1**2+XL1**2) 

GO  TO  102 

41  IF  ( XLS ( I  I  -XL )  43.43.42 

42  IM  =  i-1 
GO  TO  50 

43  DAS (II  =  ASD 

C0S(  l  »  0.0 

102  CONTINUE 

50  VOL  *  PI*R]*#2#(XLl/3.0+XL2) 

BM  =  PRHO*VOL 
400  RETURN 
STOP 
END 

COMPILE  RUN  FORTRAN 

SUBROUTINE  VELO  < XL .DXL.COS .DAS.BM. 1M.VI »DT *H»S1G»TRH0*FS»VS* 

its.tt.x.xp.vv.ff.xl: .ri.xmi 

DIMENSION  XL S( 200) .COS (200  I .DAS  I  200) 

T  *  0.0 
TT=T/TS 
V  «  VIM2.0 
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VV  «  V/(VS*12.0) 

FF-0.0 
X  -  0.0 
XLS(l)  •  0.0 
REWIND  2 

DO  410  K  •  1.1000 
IFU-H-XL1H7.25.25 

17  Fl  -  0.0 
F2  -  0.0 

DO  415  J  ■  1 . t  M 
IF  (V)  24.24.18 

18  IF ( X-H) 19.19.20 

19  I  ■  J 
XP  ■  X 
GO  TO  21 

20  IL  ■  (X-HJ/DXL 
I  ■  J+IL 

XP  »  H 

IF< I-1MJ21. 21.23 

21  IF  (XLS(I )-X)  £2.22.23 

C  COMPRESSIVE  RESISTANCE  OF  PLATE  MATERIAL 

22  DF1  -  SIG*COS< I )*DAS( I ) 

C  INERTIAL  RESISTANCE  OF  PLATE  MATERIAL 

DF2  *  0.5*TRH0*<C0$(  I  ))**$*(  VI**2*DAS<I  1 
Fl  *  F1+DF1 
F2  «=  F2+DF2 
F  «  F1+F2 

XLSU  +  1)  ■  XLS(I)+DXL 
415  CONTINUE 

23  WRITE  TAPE  2  .TT .XP .X .VV .FF 
DX=V*DT 

C  IMPULSE  MOMENTUM  LAW 

DV  «-F  *DX/IBM*V) 

V  *  V+DV 
X  »  X+DX 
T  »  T+DT 
FF  »  F/FS 
TT  »  T/T5 
VVaV/(VS*12.0) 

410  CONTINUE 

24  VaO.O 

VV  *  V/(VS*12.0) 

25  VR  »  V/12.0 

WRITE  TAPE  2.TT.XP.X.VV.FF 

XPM  ■  XP 

XM  *  X 

TM  »  T 

TTM  ■  TT 

KM-K 

PRINT  950.  VR.XPM.XM.TM.TTM 

950  FORMAT! 26H0R. VELOCITY  OF  BULLET  VR-  .E13.4.7H  FT/SEC  U 

1  26H  MAX.PtNETRAT ION  XPM  -  .E13.4.7H  INCHES  // 

2  26H  TOTAL  TRAVEL  XM  ■  .E13.4.7H  INCHES  // 

3  26H  TOTAL  TIME  FOR  XM  TM  ■  .E13.4.4H  SEC  // 

4  26H  SCALED  TIME  FOR  TM  TTM  -  .F10.3  / ) 

RETURN 

STOP 

END 
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